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Chapter  1 

INTRODUCTION 

In  the  computer  analysis  involving  wave  propagation  within  an  infinite  medium, 
a  finite  portion  of  the  medium  must  be  isolated  so  that  it  may  be  numerically  mod¬ 
elled  through  the  use  of  finite  elements.  This  creates  boundaries  at  the  perimeter  of 
the  medium  that  can  reflect  the  incident  stress  waves.  These  reflections  can  propa¬ 
gate  back  into  the  interior  of  the  medium  which  is  contradictory  to  the  actual  physical 
case.  Therefore,  it  is  appropriate  to  eliminate  or  significantly  delay  the  reflection  of 
these  waves  during  a  computer  analysis.  This  can  be  accomplished  in  either  of  two 
ways.  First,  the  dimensions  of  the  discretized  region,  or  mesh,  can  be  extended  with 
boundaries  placed  at  a  sufficient  distance  away  from  the  area  under  consideration  such 
that  the  required  response  is  obtained  before  the  reflections  return  and  contaminate 
the  solution.  This  necessitates  the  use  of  a  large  number  of  elements  which  engenders 
large  storage  requirements  and  computational  time  for  the  computer  and  may  be  un¬ 
satisfactory  to  the  user.  The  second,  which  is  the  focus  of  this  report,  is  the  use  of  a 
nonreflecting  boundary. 

A  nonreflecting  boundary  is  also  known  as  a  transmitting,  absorbing,  or  silent 
boundary.  The  purpose  of  the  nonreflecting  boundary  is  to  effectively  absorb  the 
incident  stress  waves  and,  thus,  simulate  an  infinite  space.  There  are  severed  advantages 
in  using  the  nonreflecting  boundary  in  comparison  to  the  use  of  an  extended  or  large 
mesh  with  a  fixed/free  boundary  condition.  First,  the  size  of  the  mesh  can  be  reduced 
which  decreases  the  required  amount  of  computer  storage  and  computation  time  for 
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analyses.  Along  with  the  smaller  mesh,  there  will  be  fewer  input  instructions  that  the 
user  will  need  to  prepare  or  modify.  Finally,  the  solution  can  be  performed  for  any 
length  of  time  without  the  interruption  of  spurious  reflections. 

There  are  several  beneficial  characteristics  a  nonreflecting  boundary  should  pos¬ 
sess.  In  order  for  the  nonreflecting  boundary  to  serve  its  main  purpose,  it  should  not 
require  an  excessive  amount  of  computer  storage  or  central  processing  unit  (CPU) 
time  for  it  to  function.  Likewise,  the  necessary  input  instructions  for  the  nonreflect¬ 
ing  boundary  should  be  minimal  to  save  on  user  time.  It  is  also  advantageous  for 
the  boundary  to  be  frequency  independent.  A  frequency  independent  boundary  dlows 
nonlinear  transient  analyses  to  be  easily  performed  in  the  time  domain.  In  addition,  a 
frequency-  dependent  boundary  cannot  be  used  in  an  explicit  time  solution  because  the 
natural  frequencies  of  a  system  cannot  be  calculated.  Versatility  is  another  attractive 
property  of  the  nonreflecting  boundary.  The  nonreflecting  boundary  should  be  able  to 
absorb  both  elastic  and  inelastic  stress  waves  and  be  effective  in  an  anisotropic  medium. 
These  become  important  factors  when  the  medium  is  stratified  and  highly  nonlinear, 
such  as  soil.  Although  not  a  stringent  characteristic,  the  nonreflecting  boundary  should 
be  relatively  simple  to  incorporate  into  a  computer  code. 

The  Air  Force  Weapons  Laboratory  (AFWL)  has  a  program  called  SAMSON2 
which  is  used  for  analyses  of  large  displacement,  large  strain,  nonlinear  problems  with 
the  capability  for  a  user  to  model  a  structure/media  interface.  SAMS0N2  is  the 
second  edition  of  the  SMI  (SAMSON)  code  that  was  developed  at  the  Illinois  Institute 
of  Technology  Research  Institute  (IITRI)  in  1972-73.  The  SAMS0N2  code  has  been 
subject  to  many  changes  since  its  introduction  to  the  AFWL.  There  are  intrinsic  errors 
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that  exist  in  the  code  which  are  continually  being  corrected.  Also,  the  logic  in  certain 
sections  of  SAMSON2  has  been  created  or  modified  to  increase  the  efficiency  of  the 
code. 

The  personnel  at  the  AFWL  primarily  use  SAMSON2  to  solve  for  the  dynamic 
response  of  buried  protective  structures  subject  to  blast  and/or  shock  loads.  In  order 
to  model  this  type  of  soil-structure  interaction  problem,  an  isolated  portion  of  the  sur¬ 
rounding  soil  in  addition  to  the  buried  structure  is  discretized.  Presently,  the  personnel 
at  the  AFWL  need  to  generate  a  large  mesh  with  boundaries  far  enough  away  from  the 
structure  to  allow  sufficient  time  for  the  analysis  before  the  stress  waves  return  from 
the  boundary.  A  solution  using  a  large  mesh  is  not  economical  and  requires  a  cumber¬ 
some  input  file.  To  allow  the  use  of  a  smaller  more  efficient  mesh,  it  was  desirable  to 
integrate  into  the  SAMSON2  code  a  nonrefiecting  boundary. 

The  objective  of  this  research  was  to  implement  into  the  SAMSON2  code  a 
.lonrefiecting  boundary  that  efficiently  absorbs  both  linear  and  nonlinear  stress  waves 
in  an  anisotropic  medium.  It  was  also  necessary  for  the  boundary  to  be  compatible 
with  the  SAMSON2  code.  In  other  words,  the  nonrefiecting  boundary  algorithm  must 
be  able  to  function  within  the  framework  of  the  SAMSON2  formulation. 

This  research  effort  began  with  a  review  of  the  various  nonrefiecting  boundaries 
that  exist  in  the  literature.  Of  these  boundaries,  one  needed  to  be  selected  that  was 
best  suited  for  the  SAMSON2  code.  If  none  of  the  boundaries  were  both  efficient  and 
compatible  with  SAMSON2,  a  new  nonrefiecting  boundary  had  to  be  developed.  Once 
the  decision  was  made  on  a  particular  boundary,  an  algorithm  for  the  nonrefiecting 
boundary  was  to  be  written  and  implemented  into  the  code.  The  major  portion  of 
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the  research  consisted  of  examining  several  sample  problems,  varying  in  complexity, 
that  exhibit  the  versatility  and  reliability  of  the  nonreflecting  boundary.  Though  not 
originally  part  of  the  objectives  of  the  research,  any  possible  anomalies  in  the  code  that 
were  discovered  were  to  be  noted  with  suggested  corrections.  This  report  is  a  compi¬ 
lation  of  the  work  completed  and  gives  a  summary  of  the  capabilities  and  limitations 
of  the  nonreflecting  boundary  that  weis  implemented  into  the  SAMSON2  code  . 


Chapter  2 

BACKGROUND 


An  overview  of  the  SAMSON2  code  is  necessary  in  order  to  understand  the  type 
of  nonreflecting  boundary  that  could  be  used.  The  available  nonreflecting  boundaries 
that  exist  in  the  literature  must  be  compatible  with  the  code  along  with  being  efficient 
in  the  absorption  of  the  radiant  wave  energy.  The  pertinent  characteristics  of  the 
SAMSON2  code  as  they  relate  to  a  nonreflecting  boundary  are  discussed  in  the  first 
part  of  this  chapter.  The  second  part  gives  a  brief  summary  of  the  nonreflecting 
boundaries  that  have  been  previously  developed  by  various  authors. 

2.1  The  SAMSON2  Code 

There  are  several  useful  features  within  the  SAMSON2  code.  Those  features 
that  play  an  important  role  in  the  implementation  of  a  nonreflecting  boundary  are 
presented. 

The  SAMSON2  code  is  an  efficient  two-dimensional  dynamic  finite  element  pro¬ 
gram.  The  code  utilizes  an  explicit  time  integration  scheme  in  its  formulation.  A 
central  difference  method  is  used  to  solve  for  the  displacements  and  velocities  of  each 
node  from  their  corresponding  accelerations.  The  accelerations  are  simply  computed 
from  the  equation  of  motion  at  each  node  during  a  specified  time  interval.  Therefore,  a 
system  of  equations  does  not  have  to  be  solved  simultaneously  during  each  time  step. 
This  implies  that  a  stiffness  matrix  need  not  be  generated  in  the  calculations.  The  lack 
of  a  stiffness  matrix  means  that  one  is  unable  to  calculate  the  natural  frequencies  of  a 


system  with  SAMS0N2.  Without  the  frequencies,  any  nonreflecting  boundary  that  is 
dependent  on  the  frequencies  of  the  system  cannot  be  employed. 

The  use  of  explicit  time  integration  can  create  an  instability  m  the  solution. 
The  time  step  is  restricted  to  a  maximum  length  to  maintain  stability.  This  maximum 
time  step  is  called  the  stability  limit.  The  stability  limit  is  approximately  the  time  it 
takes  for  the  fastest  stress  wave  to  travel  the  minimum  length  of  any  element  divided 
by  some  constant  factor  usually  greater  than  four. 

.Almost  all  of  the  arrays  used  in  SAMSON2  are  stored  in  a  single  large  array, 
Q.  Storage  pointers  or  indices  indicate  the  starting  location  of  each  major  array  within 
the  Q  array.  A  preprocessing  phase  of  SAMSON2  reads  the  input  parameters  while 
the  array  storage  pointers  are  saved.  For  a  particular  problem,  the  size  of  the  Q  array 
is  adjusted  before  the  actual  solution  phase  is  entered  in  order  to  reduce  the  storage 
required  by  the  parameters.  Also,  this  dynamic  storage  allocation  allows  the  size  of 
the  individual  arrays  to  be  of  any  arbitrary  length. 

The  SAMSON2  code  possesses  six  material  models  in  its  formulation.  In  addi¬ 
tion,  there  are  three  experimental  models  in  the  code  that  still  require  verification  for 
the  AFWL:  an  endochronic  model,  a  viscoplastic  model  and  a  cap  model  developed 
by  Weidlinger  Associates.  There  ate  four  elastic-plastic  material  laws  and  one  purely 
elastic  plane  strain  material  law.  The  personnel  at  the  AFWL  primarily  use  an  “en¬ 
gineering”  materird  law  in  their  soil-structure  problems.  The  “engineering”  model  is 
defined  by  a  segmented  constitutive  relationship  between  the  hydrostatic  pressure  and 
the  volumetric  strain.  Each  linear  segment  hw  an  associated  bulk  modulus  and  Pois¬ 
son’s  ratio  depending  on  whether  the  material  is  loading,  unloading,  or  reloading.  The 
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“engineering”  model  also  contains  a  yield  surface  chat  correlates  the  hydrostatic  pres¬ 
sure  to  the  square  root  of  the  second  invariant  of  the  deviatoric  stress  tensor,  yj  J',. 
If  the  current  \/~j[  value  calculated  from  the  computed  stresses  exceeds  the  failure 
surface  criterion,  the  elemental  stresses  are  adjusted  accordingly. 

.4  unique  capability  of  the  SAMSON2  code  is  a  slideline  interface  routine  that 
allows  sliding  or  separation  of  nodes  along  a  specified  line.  The  slideline  interface  is 
defined  by  master  nodes  and  slave  nodes  that  lie  along  the  interface.  The  master  nodes, 
generally,  are  nodes  on  the  coarse  side  of  the  mesh  while  the  slave  node  are  on  the  fine 
side.  The  slideline  interface  is  also  useful  in  rigidly  connecting  nodes  along  a  line  that 
separates  a  fine  mesh  from  a  coarse  mesh  without  the  need  for  transitional  elements. 
There  are  three  frictional  models  av2ulablefor  the  slideline  interface.  A  discussion  with 
Mr.  Rudeen  [12],  however,  revealed  that  the  personnel  at  the  AFWL  are  confident  in 
the  accuracy  of  only  one  of  the  frictional  laws,  the  Coulombic  sliding  friction  law.  The 
other  two,  the  impedence  law  and  the  shear  wave  transmission  law,  are  not  completely 
tested  and  verified. 

Another  unique  feature  of  the  SAMSON2  code  is  an  ability  to  calculate  a  static 
solution  using  a  technique  cadled  dynamic  relaxation.  Dynamic  relaxation  effectively 
dampens  the  lower  frequency  oscillations  that  occur  in  a  dynamic  solution.  The  lower 
frequency  modes  have  the  greatest  effect  on  the  response  of  a  system.  The  amount  of 
damping  is  controlled  by  a  mass-proportional  damping  factor  that  is  designated  in  the 
SAMSON2  input  file. 

The  SAMSON2  finite  element  library  is  divided  into  three  main  categories:  flex¬ 
ural,  continuum,  and  miscellaneous  elements.  Of  the  three  categories,  the  continuum 
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elements  are  the  ones  that  have  been  considered  in  this  study.  The  continuum  elements 
can  be  further  divided  into  triangular  and  quadrilateral  elements.  The  purpose  of  the 
triangular  element  is  to  provide  a  transition  between  the  quadrilateral  elements  and 
are  seldom  used  solely  in  an  analysis.  Therefore,  the  triangular  elements  have  not  been 
examined  extensively  in  this  study.  The  remaining  quadrilateral  elements  consist  of 
the  eight-node  isoparametric  element  (8NQ)  and  four-node  linear  element  (4NQ).  The 
personnel  at  the  AFWL  ate  not  completely  confident  in  the  results  given  by  the  8NQ. 
Also,  the  mesh  generation  capabilities  within  SAMSON2  result  in  a  more  cumbersome 
input  file  for  the  8NQ  in  comparison  to  the  4NQ.  The  4NQ,  besides  being  relatively 
easy  to  generate,  is  primarily  used  by  the  personnel  at  the  AFWL  and  is  trusted  to 
giv"  accurate  results.  Hence,  the  emphasis  in  this  study  has  been  on  the  4NQ. 

The  disadvantage  of  the  4NQ  element  is  the  tendency  for  the  element  to  exhibit 
hourglassing.  Hourglassing  is  a  mesh  instability  generated  by  an  eigenvector  displace¬ 
ment  which  results  in  no  change  in  element  forces,  that  can  result  in  large  fiuctuations 
in  the  displacements  of  the  nodes.  These  fluctuations  can,  in  turn,  make  the  element 
have  a  negative  area,  which  will  cease  execution  of  SAMSON2.  The  instability  arises 
from  one-point  quadrature  used  in  the  integration  process  within  the  element.  Fortu¬ 
nately,  an  antihourglassing  feature  is  used  to  dampen  the  effects  of  the  instability.  An 
hourglass  damping  coefficient,  tig,  is  used  to  control  the  magnitude  of  the  antihour- 
gleissing  damping  forces.  At  low  levels,  these  damping  forces  will  usually  not  affect  the 
solution  a  great  deal.  Unless  hourglassing  poses  a  problem,  fig  is  generally  set  to  zero. 
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2.2  Review  of  the  Literature 

Various  methods  have  been  developed  in  attempts  to  model  an  infinite  medium; 
all  of  which  possess  certain  limitations.  The  objective  in  each  case  is  to  avoid  little  or 
no  reflections  while  being  efficient  and  practical. 

Lysmer  and  Kuhlemeyer  [7]  developed  a  viscous  boundary  that  absorbs  the 
radiated  wave  energy  by  applying  the  following  set  of  normal  and  tangential  stresses 
to  an  otherwise  free  boundary; 

CTn  =  apc^iln  (2.1) 

r  =  bpctii,  (2.2) 

The  wave  energy  will  be  absorbed  if  the  applied  viscous  stresses  are  equal  and  opposite 
to  the  stresses  caused  by  the  incident  wave.  The  surface  stresses,  (Tn  and  t,  would 
be  proportional  to  the  dilatational  and  shear  impedences  of  the  medium,  pc^  and  pc,, 
the  normal  and  tangential  velocities  on  the  boundary,  and  U(,  and  dimensionless 
parameters  o  and  6.  The  value  for  both  of  the  dimensionless  parameters  that  gave  the 
best  overall  results  was  found  to  be  one.  The  viscous  boundary  is  a  very  good  absorber 
of  elastic  dilatational  (P)  and  elastic  shear  (S)  waves.  One  of  the  best  attributes  of  this 
boundary  is  that  it  is  independent  of  the  frequency  of  the  transmitted  wave  and  can  be 
used  within  time  domain  solutions.  Also,  it  would  be  suitable  for  both  harmonic  and 
transient  analyses.  Because  of  the  simplicity  of  the  viscous  boundary  scheme,  it  would 
be  relatively  easy  to  implement  into  a  computer  code.  The  viscous  boundary  is  local  in 
that  it  is  based  solely  on  the  characteristic  impedence  and  the  normal  and  tangential 
velocities  near  the  boundary.  This  local  attribute  becomes  important  in  the  analyses  of 


soil-  structure  interaction  problems.  The  interior  mesh  can  be  highly  nonlinear,  which 
is  a  characteristic  of  soil,  while  the  region  away  from  the  disturbance  can  remain  linear. 
The  authors  limited  the  use  of  this  viscous  boundary  to  a  plane  isotropic  medium. 

A  viscous  boundary  has  been  implemented  into  the  explicit  computer  program 
FLEX  by  Weidlinger  Associates  [18,19].  FLEX  has  been  used  for  several  applications 
since  its  introduction.  The  applications  include  seismic  wave  propagation,  electromag¬ 
netic  wave  propagation,  earthquake  engineering,  and  weapons  effects.  The  formulation 
in  FLEX  contains  a  variety  of  features  that  are  very  similar  to  those  used  by  SAM- 
SON2. 

It  is  a  popularly  held  belief  that  one  of  the  disadvantages  of  the  viscous  bound¬ 
ary  is  its  inability  to  absorb  surface  waves;  namely,  Rayleigh  waves  and  Love  waves. 
Several  authors  have  had  varied  results  when  applying  the  viscous  boundary  to  surface 
waves.  Lysmer  and  Kuhlemeyer  [7]  attempted  to  include  the  absorption  of  Rayleigh 
waves  in  their  model.  They  suggested  that  the  strength  of  the  viscous  surface  stresses 
should  also  depend  on  the  ratio  of  the  depth  from  the  free  surface  to  the  wavelength 
of  the  Rayleigh  wave.  This,  however,  makes  the  viscous  boundary  dependent  on  fre¬ 
quency  and  unsuitable  for  time  domain  solutions.  The  severity  to  which  surface  waves 
affect  the  solution  is  generally  unknown.  As  the  dilatationsd  waves  and  shear  waves 
propagate  outward,  they  form  surface  waves  as  they  strike  the  boundary  at  large  angles 
of  incidence.  This  may  not  necessarily  be  undesirable.  If  surface  waves  are  generated, 
they  will  travel  relatively  slowly  along  the  boundary  without  reffecting  back  into  the 


interior  of  the  medium. 


If  the  directions  of  the  incident  stress  waves  were  known,  the  viscous  boundary 
would  completely  absorb  the  waves.  Unfortunately,  only  in  a  one-dimensional  wave 
propagation  problem  is  the  angle  of  incidence  known.  In  two-  or  three-dimensional 
problems  the  angle  of  incidence  varies  between  0  degrees  (wave  impinging  normal  to 
the  surface)  and  90  degrees  (wave  travelling  tangentially  to  the  surface).  The  efficiency 
of  the  viscous  boundary  in  absorbing  the  stress  waves  is  dependent  upon  the  angle  of 
incidence  of  the  incoming  stress  waves.  At  incident  angles  greater  than  about  60  degrees 
the  efficiency  of  the  viscous  boundary  begins  to  noticeably  decrease.  Usually,  the  -  ave 
energy  that  is  not  absorbed  when  the  angle  of  incidence  is  large  is  of  little  concern 
for  the  following  reasons.  First,  one  can  design  a  mesh  such  that  the  nonreflecting 
boundary  is  oriented  normal  to  the  wave  source.  Second,  those  waves  that  do  impinge 
at  large  angles  of  incidence,  usually  strike  one  or  more  edges  of  the  viscous  boundary 
at  lower  angles  of  incidence,  where  they  can  be  efficiently  absorbed,  before  radiating 
back  toward  the  source.  Also,  as  previously  mentioned,  some  of  the  wave  energy  that 
is  not  absorbed  is  transformed  into  surface  waves.  The  surface  waves  travel  along  the 
boundary  and,  therefore,  do  not  affect  the  interior  of  the  mesh. 

The  efficiency  of  the  viscous  boundary  is  also  dependent  upon  Poisson’s  ratio  of 
the  medium.  Cohen  and  Jennings  [3]  noted  that,  as  Poisson’s  ratio  tends  to  its  maxi¬ 
mum  value  of  0.5,  the  absorption  characteristics  of  the  viscous  boundary  are  generally 
less  efficient.  This  is  particularly  noticeable  when  an  incident  S-wave  reflects  a  P-wave. 
However,  for  normal  values  of  Poisson’s  ratio,  the  effect  on  the  efficiency  is  negligible. 

White,  Valliappan,  and  Lee  [20]  proposed  to  refine  the  Lysmer  and  Kuhlemeyer 
viscous  boundary  by  making  it  more  efficient  and  versatile.  An  iterative  procedure 


was  developed  to  optimize  the  absorption  qualities  of  the  boundary.  The  stresses  that 
are  not  initially  absorbed  by  the  viscous  boundary  are  minimized  with  respect  to  the 
initial  damping  matrix.  A  new  damping  matrix  is  then  calculated,  and  the  procedure 
is  repeated.  For  an  isotropic  medium,  the  optimized  values  of  the  dimensionless  pa¬ 
rameters,  a  and  6,  were  found  as  functions  of  Poisson’s  ratio.  Unfortunately,  the  use 
of  the  optimized  parameters  resulted  in  a  very  small  increase  in  the  efficiency  of  the 
boundary.  The  authors  also  stated  that  an  anisotropic  medium  could  be  modelled  by 
this  iterative  procedure.  However,  they  have  not  concluded  how  efficient  the  new  “uni¬ 
fied”  viscous  boundary  actually  is  for  an  anisotropic  medium.  In  addition,  the  authors 
performed  a  study  of  an  axisymmetric  medium.  They  discovered  that  the  efficiency 
of  a  viscous  boundary  parallel  to  the  axis  of  symmetry  is  dependent  upon  the  ratio  of 
the  distance  away  from  the  axis  of  symmetry  (R)  to  the  wavelength  of  the  stress  wave 
(A).  For  R/A  values  greater  than  0.5,  the  axisymmetric  coefficients  are  virtually  the 
same  ew  for  the  plane  strain  ceise.  For  a  viscous  boundary  perpendicular  to  the  etxis  of 
symmetry,  there  was  a  negligible  decrease  in  the  efficiency  if  plane  strain  coefficients 
were  used  instead  of  the  axisymmetric  coefficients  regardless  of  the  value  of  R/  A. 

Lysmer  and  Waas  [8]  developed  an  excellent  nonreflecting  boundary  for  the 
transmission  of  shear  waves  within  a  layered  medium.  This  boundary  has  been  suc¬ 
cessfully  used  in  the  computer  program  FLUSH.  Although  the  authors’  work  involved 
only  Love  waves,  the  application  of  this  boundary  into  FLUSH  has  been  extended  to 
also  transmit  Rayleigh  waves.  The  authors’  theory  is  to  incorporate  the  stiffness  of 
the  outer  region  into  the  equation  of  motion  at  the  boundary.  A  “boundary  stiffness” 
matrix  is  calculated  from  the  wave  numbers  of  the  incident  stress  waves.  The  boundary 
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stresses  that  are  applied  to  transmit  the  incident  waves  are  then  proportional  to  the 
displacements  of  the  nodes  on  the  boundary.  Unfortunately,  this  boundary  has  several 
limitations.  The  primary  limitation  is  that  the  boundary  cannot  transmit  dilatational 
waves.  Because  the  wavelength  of  the  impinging  wave  must  be  known  to  calculate 
the  boundary  stiffness  matrix,  the  boundary  is  frequency  dependent.  Therefore,  the 
boundary  is  not  applicable  to  time  domain  solutions  and  difficult  to  use  in  transient 
analyses.  Another  restriction  of  the  boundary  is  that  the  interior  region  of  the  mesh 
must  remain  linear.  Finally,  the  nonreflecting  boundary  can  only  be  applied  to  the 
lateral  sides  of  the  mesh  while  the  base  must  be  fixed. 

In  order  to  eliminate  all  reflections  from  the  boundary,  the  superposition  of 
solutions  that  satisfy  the  Dirichlet  and  Neumann  boundary  conditions  was  suggested 
by  Smith  [16].  In  two  dimensions,  the  superposition  boundary  is  analagous  to  a  pair 
of  regions  in  which  boundaries  of  one  region  are  fixed  in  the  normal  direction  and 
free  in  the  tangential  while  the  boundaries  for  the  other  region  are  the  opposite  to 
the  first  region.  The  addition  of  the  two  solutions  results  in  the  cancellation  of  all 
reflections,  including  both  body  waves  and  surface  waves.  This  method  restricts  the 
entire  mesh  to  remain  elastic  so  that  the  principle  of  superposition  can  be  applied.  The 
superposition  boundary  becomes  exceedingly  impractical  for  almost  ail  cases.  Each 
possible  combination  of  reflections  near  the  boundary  surface  must  be  accounted  for 
by  two  solutions.  In  other  words,  for  n  nonreflective  surfaces,  the  complete  analysis 
requires  2**  solutions.  The  use  of  the  superposition  boundary  becomes  prohibitive  even 
for  a  region  that  has  only  one  or  two  nonreflecting  surfaces. 
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Cund^dl,  et  al,  [5]  and  Kunar  and  Rodriguez-Ovejero  [6j  developed  an  idea  to 
eliiTiinate  the  problem  of  the  large  number  of  solutions  required  by  Smith's  superpo¬ 
sition  boundary.  The  objective  of  their  formulation  was  to  have  only  a  small  area  of 
the  mesh  require  a  superimposed  solution.  Two  narrow  boundary  regions  consisting 
of  a  depth  of  3  or  4  elements  would  be  connected  to  the  interior  mesh.  Instead  of 
fixed-free  boundary  conditions,  the  surfaces  of  the  boundary  regions  would  have  con¬ 
stant  velocity-constant  stress  similar  to  the  scheme  used  by  Smith.  At  every  third  or 
fourth  time  step  the  value  for  each  variable  in  the  two  boundary  regions  is  updated 
by  averaging  the  corresponding  vcilues  in  both  boundary  regions.  The  stresses  and 
velocities  on  the  boundary  surface  are  kept  constant  during  the  iaterim  time  steps. 
The  distribution  of  the  stresses  and  velocities  near  the  boundary  nearly  behaves  « 
if  the  attached  boundary  region  extended  to  infinity.  The  major  advantage  of  this 
boundary  is  the  drastic  reduction  in  the  need  for  multiple  solutions.  Furthermore,  only 
the  boundary  regions  need  to  remain  linear.  The  refined  superposition  boundary  will 
cancel  both  body  and  surface  waves  regardless  of  the  angles  of  incidence.  The  superpo¬ 
sition  boundary  has  been  successfully  used  in  the  finite  difference  program  DAMSEL. 
Perhaps  the  only  disadvantages  to  this  boundary  are  that  it  effectively  adds  a  layer 
of  6  or  8  elements  on  the  boundary  of  the  mesh  whose  number  may  be  significant  in 
comparison  to  the  number  of  elements  in  the  interior  mesh.  Also,  the  variables  within 
the  boundary  regions  must  be  saved  throughout  the  calculations. 

Robinson  [10]  separated  the  dilatational  waves  and  shear  waves  into  their  two 
respective  potentials  near  the  vicinity  of  the  transmitting  boundary.  The  author  as¬ 
sumed  that  the  body  waves  travel  in  both  a  plane  and  a  cylindrical  pattern.  For  the 


plane-wave  approximation,  the  transmitting  boundary  gave  adequate  results  only  at 
points  away  from  the  boundary.  However,  when  the  cylindrical  wave  approximation 
was  used,  the  boundary  yielded  a  very  accurate  solution  for  the  entire  region.  The 
least  accurate  results  in  both  cases  appear  at  the  junction  of  the  axis  of  symmetry 
and  the  transmitting  boundary.  There  is  no  restriction  on  the  linearity  of  the  inte¬ 
rior  mesh;  only  the  region  near  the  transmitting  boundary  needs  to  remain  linear. 
The  cylindrical-wave  approximation  becomes  more  accurate  with  an  increeise  in  the 
curvature  of  the  wave.  The  need  for  a  large  curvature  (or  large  radius)  of  the  stress 
wave  at  the  boundary  implies  that  the  transmitting  boundary  must  not  be  placed  too 
close  to  the  wave  source  in  order  to  obtain  an  accurate  solution.  Therefore,  the  size 
of  the  mesh,  in  some  cases,  may  still  be  too  large  to  allow  the  practical  use  of  this 
nonreflecting  boundary. 

Cohen  and  Jennings  [3]  made  improvements  to  a  paraxial  boundary  developed 
from  previous  authors  which  they  called  the  extended-parrudal  boundary.  A  paraxial 
medium,  which  transmits  waves  in  one  direction  only,  is  attached  to  the  interior  mesh. 
There  were  some  obstacles  to  overcome  before  the  extended-parajcial  boundary  could 
be  effectively  implemented.  An  interface  element  linking  the  paraxial  boundary  and 
the  interior  elastic  mesh  has  to  be  used  for  a  smooth  wave  transmission  between  the  two 
regions.  Because  the  paraxial  boundary  creates  a  nonsymmetrical  boundary  damping 
matrix,  spurious  oscillations  will  result  in  the  solution.  A  numerical  dissipative  proce¬ 
dure  is  used  to  eliminate  the  undesired  oscillations.  Instabilities  were  observed  when 
Poisson’s  ratio  of  the  medium  was  greater  than  1/3.  This  induced  a  negative  term 
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in  the  stiffness  matrix.  For  Poisson’s  ratio  greater  than  1  '3,  the  authors  equated  this 
relatively  insignificant  negative  term  to  zero  which  eliminated  the  instability. 

A  thorough  study  of  the  extended-paraxial  boundary  and  the  viscous  bound¬ 
ary  was  made  by  Cohen  and  Jennings  to  compare  the  various  characteristics  of  each 
boundary.  It  was  noted  that  a  viscous  boundary  is  simply  the  first  order  paraxial 
boundary.  The  same  advantages  exist  in  both  boundaries  except  that  the  extended- 
paraxial  boundary  is  theoretically  more  accurate  than  the  viscous  boundary.  However, 
in  actual  numerical  results,  this  superiority  is  not  significant.  The  authors  showed  that 
the  extended-paraxial  boundary  is  more  efficient  in  transmitting  Rayleigh  waves  than 
the  viscous  boundary.  The  viscous  boundary,  however,  still  performed  adequately  in 
the  absorption  of  the  Rayleigh  waves.  Similar  to  the  work  done  by  White,  et  al,  Cohen 
and  Jennings  found  that  the  viscous  boundary  parallel  to  the  axis  of  symmetry  does 
become  less  efficient  for  lower  v^dues  of  R/A.  Both  the  viscous  and  extended- paraxial 
boundaries  did  not  absorb  high  frequency  waves  very  effectively.  The  latter  actually 
amplified  the  noise.  These  oscillations  are  ewily  eliminated  by  using  some  numerical 
damping  that  filters  out  only  high  frequency  waves.  A  characteristic  not  pointed  out 
by  the  previous  authors  concerning  the  viscous  boundary  is  the  insensitivity  of  the 
absorption  efficiency  to  the  parameters  a  and  h.  This  implies  that  the  viscous  bound¬ 
ary  may  be  able  to  absorb  stress  waves  of  varying  velocities  including  slower  travelling 
nonlinear  waves.  In  the  numerical  application  of  the  two  nonreflecting  boundaries, 
the  authors  noted  that  neither  boundary  had  an  adverse  effect  on  the  stability  of  the 
solution.  Finally,  the  extended-paraxial  boundary  is  slightly  more  costly  to  implement 
than  the  viscous  boundary. 
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Of  the  several  nonreflecting  boundaries  presented,  only  a  few  are  not  suitable  for 
the  SAMSON2  code.  These  include  the  boundaries  that  are  frequency  dependent,  not 
compatible  with  the  code,  or  do  not  allow  the  original  size  of  the  mesh  to  be  reduced 
by  an  appreciable  degree.  The  nonreflecting  boundary  that  would  be  implemented  into 
the  SAMSON2  code  would  be  chosen  based  on  the  boundary’s  efficiency  of  transmission 
versus  the  storage  requirements  and  computation  time  for  the  computer  and  the  ease 
of  implementation  into  the  SAMSON2  code. 


IS 


Chapter  3 

THE  IMPLEMENTATION  OF  THE  NONREFLECTING  BOUNDARY 


The  viscous  boundary  scheme  developed  by  Lysmer  and  Kuhlemeyer  [2]  was 
initially  chosen  to  be  implemented  into  the  SAMSON2  code.  Of  the  various  nonre¬ 
flecting  boundaries  that  were  investigated,  the  viscous  boundary  appeared  to  absorb 
the  wave  energy  effectively  and  seemed  to  be  particularly  well  suited  for  the  code.  The 
improvements  made  by  White,  et  al,  [6]  had  such  a  negligible  effect  on  the  efficiency 
of  the  viscous  boundary  that  they  were  ignored  in  this  study.  The  only  other  nonre¬ 
flecting  boundary  considered  in  this  study  w&s  the  extended-paraxial  boundary.  The 
extended-paraxial  boundary  is  slightly  more  costly  to  use  than  the  viscous  boundary 
and  more  difficult  to  implement.  For  these  two  reasons,  the  viscous  boundary  was  the 
preliminary  choice.  If  the  viscous  boundary  did  not  perform  to  an  acceptable  level,  the 
extended- paraxial  boundary  would  then  be  implemented. 


3.1  The  Theory  of  the  Viscous  Boundary 

The  concept  of  the  viscous  boundary  can  be  readily  understood  by  examining 
the  case  of  one-dimensional  wave  propagation  in  a  bar.  The  governing  differential 
equation  for  plane  strain  is 


d^u  _  d^u 


(3.1) 


where,  u  =  horizontal  displacement  variable 

p  =  mass  density 
E  =  Young’s  modulus 
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The  solution  of  Equation  3.1  is, 

u  =  f(ct  ~  x)  +  g{ci  +  x) 


(3,2) 


where, 


(3.3) 


The  arbitrary  functions,  /  and  g,  represent  two  waves:  one  propagating  in  the  positive 
direction  and  one  propagating  in  the  negative  direction.  The  exact  waveshapes  are 
determined  from  the  boundary  conditions  and  the  initial  conditions  of  the  problem. 
Assume  that  x  is  measured  from  the  right  end  of  a  bar  and  that  the  right  end  is  free  to 
translate  in  the  i-direction.  For  a  wave  travelling  in  the  negative  i-direction  (to  the 
right),  the  displacement  function  of  the  incident  wave  is  given  by  the  function  g.  The 
displacement  function  for  the  wave  reflected  from  the  free  end  is  given  by  the  function 
/.  Therefore,  the  stress  due  to  the  incident  wave  is 

^  ^  ~  +  *)  (3.4) 

and  the  stress  due  to  the  reflected  wave  is 


<Tref=E^  =  -Ef'{ct~x)  (3.5) 

Assume  a  stress,  <r„,  is  applied  to  the  free  right  end  of  the  bar.  The  summation  of  the 
stresses  at  the  right  end  must  equal  zero. 


(Tin  +  cr,*/  +  o-v  =  0  (3.6) 

In  order  for  the  reflected  wave  stress  in  Equation  3.6  to  be  zero, 

<^in  —  ~^v 


(3.7) 
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Because  the  applied  stress  a,  is  proportional  to  the  velocity  at  the  end  of  the  bar,  the 
stress  is  viscous  and  will  dissipate  the  energy  of  the  incident  stress  wave.  The  constant 
term,  pc,  is  called  the  characteristic  impedence  of  the  material.  The  negative  sign 
implies  that  the  viscous  stress  must  be  applied  opposite  in  direction  to  the  velocity  at 
the  boundary. 

The  simplicity  of  the  one-dimensional  wave  propagation  problem,  unfortunately, 
does  not  exist  in  problems  of  two  or  more  dimensions.  Figure  3.1  shows  an  incident 
dilatational  wave  (P-wave)  impinging  on  a  boundary  at  an  angle  a.  The  P-wave  will 
reflect  as  another  P-wave  at  the  same  angle  a  and  as  a  S-wave  at  an  angle  of  (3.  The 
similar  case  of  an  incident  S-  wave  is  shown  in  Figure  3.2. 

The  viscous  boundary  scheme  proposed  by  Lysmer  and  Kuhlemeyer  equated  the 
normal  and  tangential  stresses  of  Equations  2.1  and  2.2  to  the  stresses  at  the  boundary 
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S-wave 


Figure  3.1  Incident  P-wave 


Figure  3.2  Incident  S-wave 
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caused  by  the  incidence  of  the  P-wave  and  S-wave.  The  viscous  stresses  will  perfectly 
absorb  the  incident  stress  waves  if  they  impinge  normal  to<he  boundary.  The  efficiency 
of  the  viscous  boundary,  however,  decreases  as  the  angle  of  incidence  of  both  waves 
increases  from  zero. 


3.2  The  Development  of  the  Equations  for  the  Nonreflecting  Boundary 

The  wave  speeds  of  the  P-wave  and  S-wave  and  the  mass  density  of  the  medium 
dictate  the  value  of  the  dilatational  and  shear  impedences.  The  relationship  between 
the  two  impedences  depends  on  the  value  of  Poisson’s  ratio.  The  dilatational  wave 
speed  is  given  by 

c,  =  (3.12) 


where,  M  =  constrained  modulus 

p  =  mass  density 
The  shear  wave  speed  is  given  by 


(3.13) 


where,  G  =  shear  modulus 

The  modified  formulation  of  the  viscous  boundary  in  the  SAMSON2  code  is  based 
on  Young’s  modulus.  From  the  following  two  equations  that  relate  the  constrruned 
modulus  and  shear  modulus  to  Young’s  modulus  by  the  Poisson's  ratio,  i/, 


(3.14) 


G  = 


2(1  +  0 


(3.15) 
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From  the  impedences  of  the  elements  near  the  boundary,  the  viscous  forces  at 
the  boundary  can  be  calculated.  Figure  3.3  shows  a  typical  element  (in  this  case, 
a  four-node  quadrilateral  element)  adjacent  to  a  nonreflecting  boundary.  The  nodes 
i  and  i  +  1  are  nodes  that  define  a  segment  along  the  nonreflecting  boundary.  The 
horizontal  and  vertical  velocities  of  the  nodes  ate  represented  as  u  and  v,  respectively. 
He  is  the  horizontal  projection  and  Vf  is  the  vertical  projection  of  the  line  connecting 
nodes  t  and  t  -t-  1  to  the  global  x  and  y  axes. 

The  normal  and  tangential  stresses  applied  on  the  boundary  can  be  transformed 
into  their  respective  horizontal  and  vertical  components.  The  horizontal  and  vertical 
viscous  stresses  applied  at  the  nodes,  therefore,  are  both  dependent  on  the  contributing 
effects  of  the  impinging  S-waves  and  P-waves.  For  simplicity,  the  case  of  an  incident 
plane  P-wave  is  examined.  The  horizontal  component  of  the  normal  viscous  stress  is 


<Tx  =  -pCdUx 


(3.20) 
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V,rl 

U|  I 


where,  it,  is  the  horizontal  component  of  the  normal  velocity. 


it,  =  +  Nt+i{y)u,+i 

(3.21) 

where  lV,(y)  and  N^^i{y)  are  the  linear  interpolation  functions  of  the  horizontal  com¬ 
ponent  of  the  normal  velocity  whose  origin  is  at  node  t 

(3.22) 

'V.ti(y)  =  ^ 

(3.23) 

The  horizontal  forces  at  the  nodes  can  be  calculated  using  the  linear  interpolation 
functions  of  Equations  3.22  and  3.23 

Fx,=-  Nx{y)<7ztdy 

Jo 


(3.24) 


where 


Fz.^i  =  -  /  -Vi-ri(!/)(7iiciy  (3.25) 

Jo 


,  t  =  thickness  of  the  element  normal  to  the  plane  of  the  element 

Substituting  Equations  3.20  and  3.21  into  Equations  3.24  and  3.25 

-V. 


f'x,  =  -  f  PCjN^(y)[N^(y)u^  .V,  +  ,(y)^^  +  l] 

Jo 

=  -  /  pCd-^i+iiy)[Nx{y)ut  +  iV,+i(y)ut+i] 
Jo 


(3.26) 

(3.27) 


Integrating  and  putting  into  matrix  form, 


(3.28) 


By  a  similar  argument,  the  horizontal  nodal  forces  caused  by  an  incident  plane  S- 
wave  can  be  found.  In  this  case,  the  linear  interpolation  functions  for  the  horizontal 
component  of  the  tangenti£d  velocity  must  now  be  used. 


Ar.(x)  =  1  -  I; 

=  d; 

Resulting  in  the  equation 


pC/tHf 


2/6  1/6 
1/6  2/6 


(3.29) 

(3.30) 


(3.31) 


Combining  the  total  horizontal  forces  from  the  incident  P-wave  and  the  incident  S- 
wave. 


{  J  [2(6  (3.32) 
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Expanding  this  procedure  further,  the  total  vertical  forces  at  the  two  nodes  are 


-pcd^Ht 


’2/6 

1/6 


(3.33) 


Equations  3.32  and  3.33  are  the  final  versions  of  the  equations  that  are  implemented 
into  the  SAMSON2  code.  The  horizontal  viscous  force  vector  including  Fx,  and 
is  represented  by  FXl  and  FX2  in  the  code.  Likewise,  the  vertical  viscous  force  vector 
is  represented  by  FYl  and  FY2.  Also,  the  vertical  distance  and  the  horizontal  distance, 
Vf  and  Ht,  are  expressed  eis  CC  and  SS  respectively.  The  variable  name  THICK  is  the 
thickness  t. 

The  preceding  equations  assumed  planar  stress  waves.  For  the  axisymmetric 
caise,  all  of  the  forces  are  multiplied  by  the  distance  the  midpoint  of  the  segment 
joining  nodes  i  and  t  +  1  is  from  the  axis  of  symmetry. 


3.3  The  Formulation  of  the  Nonreflecting  Boundary  in  the  SAMSON2  Code 
The  riscous  boundary,  as  originally  developed,  is  limited  to  an  isotropic  plane 
medium.  Also,  the  region  near  the  boundary  must  remain  elastic.  The  features  within 
the  SAMSON2  code  provided  the  means  to  extend  the  capability  of  the  viscous  bound¬ 
ary.  The  wave  energy  propagating  in  an  anisotropic,  axisymmetric  medium  could  also 
be  absorbed  by  using  the  modified  formulation  within  the  SAMSON2  code.  The  effec¬ 
tiveness  of  the  viscous  boundary  in  absorbing  inelastic  waves,  however,  was  unknown. 
Cohen  and  Jennings  [3]  stated  that  the  viscous  boundary  could  be  applied  for  inelastic 
problems  due  to  the  insensitivity  of  the  viscous  stresses  to  the  dimensionless  parame¬ 
ters  a  and  b.  However,  they  did  not  solve  any  problems  that  had  nonlinear  behavior. 


Several  problems  were  run  in  this  study  to  test  the  effectiveness  of  the  viscous  boundary 
in  an  inelastic  medium.  The  results  will  be  discussed  in  Chapter  4. 

The  principle  material  law  used  by  the  AFWL  is  the  “engineering”  material 
model.  During  each  time  step,  the  “engineering”  material  model  subroutine  updates 
the  stiffness  within  each  element  (i.e.,  a  nonlinear  material  law).  This  implies  th’*  the 
characteristic  impedence  for  each  element  can  also  be  updated  at  each  time  step.  For 
an  element  located  on  the  boundary,  the  strengths  of  the  viscous  stresses  applied  at  its 
boundary  are  dictated  by  its  current  impedences.  Consequently,  the  viscous  boundary 
should  be  able  to  absorb  the  radiating  wave  energy  in  an  inelastic,  anisotropic  medium. 
Unfortunately,  when  the  logic  was  programmed  into  the  SAMSON2  code  and  several 
test  runs  were  made,  the  results  were  less  accurate  than  those  that  were  based  on  a 
purely  el^lstic  viscous  boundary. 

The  strategy  for  implementing  the  nonreffecting  boundary  into  SAMSON2  was 
to  create  as  few  additional  executable  lines  of  code  as  possible.  Fewer  lines  would  limit 
the  amount  of  CPU  time  required  for  the  operation  of  the  nonreffecting  boundary. 
Thus,  there  was  no  incorporation  of  an  additional  subroutine  in  the  code.  The  existing 
subroutines  were  to  be  modiffed  to  provide  the  means  to  implement  a  nonreffecting 
boundary. 

Within  the  subroutine  FREEF2,  the  external  forces  of  the  system  are  calculated 
and  located  in  the  array  FORCD.  FREEF2  mainly  uses  input  information  from  the 
load  line  data  cards  to  calculate  these  external  forces.  Among  the  arguments  passed 
to  FREEF2  are  the  velocity  array,  V,  and  the  load  line  information  arrays:  IVOL, 
KPRES  and  PT.  IVOL  contains  the  load  line  master  control  data,  KPRES  lists  the 
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node  numbers  along  the  load  line,  and  PT  contains  either  time  function  pairs  or  the 
initial  velocity  data.  The  nonreflecting  boundary  can  be  interpreted  as  externally 
applied  forces  that  depend  on  the  velocity  of  the  nodes  on  the  boundary  and  the  local 
impcdences  of  the  medium  near  the  boundary.  The  nonreflecting  boundary,  therefore, 
can  be  treated  as  a  load  line.  Thus,  the  main  formulation  of  the  nonreflecting  boundary 
was  located  in  the  subroutine  FREEF2.  The  velocity  array  is  already  available  to  the 
subroutine  in  the  original  SAMSON2  code  (passed  to  FREEF2  as  array  V).  Only  an 
array  of  the  impedences  of  the  material  near  the  boundary  needed  to  be  added  to  the 
list  of  arrays  that  are  passed  to  FREEF2. 

The  subroutine  READPV  is  used  to  read  the  load  line  input  information.  The 
modifications  made  to  READPV  by  the  authors  would  entail  two  objectives.  First, 
the  subroutine  must  allow  the  input  instructions  for  the  nonreflecting  boundary  to  be 
read.  Second,  the  impedences  of  the  medium  near  the  boundary  would  be  calculated 
and  stored  in  an  array.  The  new  modified  version  of  READPV  is  listed  in  Appendix 
C. 

The  subroutine  READPV  is  called  NPRES  times  by  the  subroutine  PREPROC, 
where  NPRES  is  the  total  number  of  load  lines  in  a  problem.  In  the  original  SAMSON2 
code,  the  arguments  passed  to  READPV  are  the  load  line  information  arrays  and  the 
variables  THICK  and  LADD.  THICK  is  the  thickness  of  a  planar  problem  pressure 
line,  and  LADD  is  the  extra  storage  requirement  for  the  load  line  information.  In 
the  modified  version,  four  additional  arrays  are  peissed  to  READPV:  E,  INDXE,  IX 
and  BIMP.  E  contains  the  element  group  data.  The  E  array  is  split  into  three  different 
“subarrays”  for  each  material  group.  The  three  subarrays  are  the  element-group  master 
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control  array,  the  element  parameter  array  and  the  material  law  parameter  array.  The 
storage  pointers  that  list  the  starting  locations  of  the  three  subarrays  for  each  element 
group  are  stored  in  array  INDXE.  The  material  law  parameter  subarray  is  the  only 
one  that  is  utilized  in  READPV.  IX  contains  the  element  connectivity  data.  BIMP  is 
a  new  array  that  stores  the  dilatational  and  shear  impedences  for  each  element  along 
the  nonreflecting  boundary.  When  the  arrays  KPRES,  PT,  and  BIMP  are  passed  to 
READPV  from  PREPROC,  they  are  all  indexed  from  the  same  storage  pointer.  The 
call  statement  passes  these  three  arrays  as  Q(L-l-6). 

At  the  beginning  of  subroutine  READPV,  the  load  line  control  data  on  card  lOA 
are  initially  read.  This  procedure  has  not  been  changed  by  the  authors.  The  next  step 
is  an  IF-  THEN-ELSE  block  that  reads  time  function  pairs  (if  IVOL(l)  equals  -2,  -1,  1, 
or  2)  or  the  initial  velocity  data  (if  IVOL(l)  equals  0)  of  card  lOB.  The  values  are  stored 
in  the  array  PT.  An  additional  ELSEIF  statement  has  been  included  in  this  block.  If 
IVOL(l)  equals  3  or  4  (the  load  line  is  a  nonreflecting  boundary),  the  number  of  time 
function  pairs  is  zero  and  the  value  is  stored  in  IVOL(5).  The  node  numbers  along 
the  load  line  are  then  read  and  stored  in  the  array  KPRES.  Because  PT  and  KPRES 
were  indexed  at  the  same  location  when  they  were  passed  to  READPV,  the  values  of 
KPRES  are  stored  beginning  at  location  2*NPT-(-l.  This  procedure  has  also  not  been 
changed  from  the  original  code.  The  remainder  of  the  subroutine  is  the  addition  of  a 
double  DO  loop  block  that  finds  the  element  numbers  along  a  nonreflecting  boundary, 
calculates  the  dilatational  and  shear  impedences  for  each  element  and  stores  the  values 


in  the  array  BIMP. 


The  first  DO  loop  is  incremented  from  the  starting  location  of  the  array  KPRES 
to  the  second  from  last  location  of  KPRES.  The  purpose  of  the  DO  loop  is  to  retrieve 
two  successive  node  numbers  along  the  boundary  at  a  time.  Each  pair  of  successive 
node  numbers  define  a  segment  on  the  boundary  and  an  edge  of  a  particular  element. 
The  second  DO  loop  then  checks  each  element  in  the  mesh  and  finds  which  element 
has  an  edge  that  corresponds  to  the  current  pair  of  node  numbers  on  the  boundary. 
.A.t  the  beginning  of  the  DO  loop,  the  type  of  each  element  is  found  (i.e.,  the  number 
of  nodes  of  an  element).  Depending  on  the  type  of  the  element,  all  of  the  edges  of  the 
element  (two  adjacent  nodes  on  the  element)  will  be  compared  to  the  current  pair  of 
node  numbers  on  the  boundary.  The  nodes  numbers  are  checked  in  a  counterclockwise 
sense.  In  other  words,  the  element  will  be  to  the  left  of  the  line  connecting  the  first 
node  to  the  second  node.  If  none  of  the  edges  of  the  element  match  the  current  segment 
on  the  boundary,  the  next  element  is  checked.  If  an  edge  of  the  element  matches  the 
current  segment,  the  second  DO  loop  is  exited.  The  dilatational  and  shear  impedences 
of  the  element  are  then  calculated  and  stored. 

The  element  group  number  of  the  element  is  used  to  find  the  storage  pointer  of 
the  material  law  parameter  array.  From  this  array,  the  density.  Young’s  modulus  and 
Poisson’s  ratio  within  the  element  are  found.  The  dilatational  and  shear  impedences 
of  the  element  are  calculated  using  Equations  3.18  and  3.19  and  stored  in  the  array 
BIMP.  Similar  to  KPRES,  BIMP  was  passed  to  READPV  indexed  at  the  same  location 
as  PT.  As  previously  mentioned,  the  first  DO  loop  begins  from  the  storage  pointer, 
2*NPT+1,  where  the  array  KPRES  is  indexed.  The  element  impedences  must  be  stored 
starting  from  the  location  immediately  after  KPRES.  Therefore,  the  storage  pointer 
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for  the  array  BIMP  is  NDNOD  greater  than  the  storage  pointer  for  KPRES,  where 
NDNOD  is  the  number  of  nodes  on  the  load  line.  An  additional  counter  is  used  to  allow 
for  the  storage  of  two  variables  for  each  pass  of  the  first  DO  loop.  The  dilatational 
impedence  is  stored  first  and  the  shear  impedence  second.  After  the  impedences  have 
been  calculated  and  stored  for  a  segment  on  the  boundary,  the  next  pair  of  node 
numbers  is  examined,  and  the  procedure  is  repeated  until  all  of  the  nodes  on  the  load 
line  have  been  accounted  for.  The  leist  step  of  the  subroutine  is  to  adjust  the  value 
of  LADD  for  the  additional  storage  requirement  for  the  array  BIMP.  For  a  series  of 
nodes  along  a  load  line,  the  number  of  elements  will  be  one  less  than  the  number  of 
nodes.  Thus,  the  variable  LADD  must  be  increased  by  NDNOD-1  for  a  nonreflecting 
boundary. 

The  subroutine  SOLVE  controls  the  main  solution  phase  of  SAMSON2.  SOLVE 
begins  by  calling  the  subroutine  FREEFD  that  calculates  the  initial  boundary  condi¬ 
tions  of  a  problem.  The  main  time  integration  loop  is  then  entered.  Within  the  loop, 
the  external  forces  of  the  problem  are  computed  by  the  subroutine  FREEF2,  the  inter¬ 
nal  forces  by  the  subroutine  FRCIN  and  the  slideline  forces  by  the  subroutine  SLIDER. 
With  these  forces,  the  displacements,  velocities  and  accelerations  of  the  nodes  are  found 
using  the  equations  of  motion  by  calling  the  subroutine  MOTION.  There  have  only 
been  two  modifications  made  by  the  authors  to  the  subroutine  SOLVE.  The  two  mod¬ 
ifications  are  to  the  call  statements  for  FREEFD  and  FREEF2.  In  the  modified  cod^, 
the  array  BIMP  is  also  passed  to  these  two  subroutines.  The  array  is  passed  as  PR(L3). 
In  the  subroutine  SOLVE,  all  of  the  load  line  arrays  are  stored  in  the  array  PR.  The 
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variable  L3  is  the  storage  pointer  for  the  array  KPRES.  The  reason  for  the  modifi¬ 
cation  to  the  call  statement  to  subroutine  FREEFD  is  that  FREEFD  later  calls  the 
subroutine  FREEF2. 

The  changes  to  subroutine  FREEFD  were  made  solely  to  allow  for  the  call  to 
FREEF2.  The  array  BIMP  has  been  added  to  the  argument  list  of  the  first  line  of 
FREEFD.  Correspondingly,  the  size  of  the  array  must  be  declared  to  be  one  in  the 
dimension  statement.  The  last  modification  is  the  addition  of  the  array  BIMP  to  the 
arguments  passed  in  the  call  statement  to  FREEF2. 

The  new  modified  subroutine  FREEF2  is  listed  in  Appendix  C.  Only  two  lines 
in  the  subroutine  have  been  changed;  all  other  modifications  have  been  only  the  addi¬ 
tion  of  new  lines.  The  argument  list  in  the  first  line  of  FREEF2  now  lists  an  additional 
array  called  BIMP.  Recall  that  BIMP  contains  the  boundary  impcdences  of  the  ele¬ 
ments  along  a  particular  nonreflecting  boundary  that  were  calculated  in  the  subroutine 
READPV.  The  second  line  modification  is  the  dimension  statement  where  the  size  of 
the  array  BIMP  is  declared  to  be  one. 

Upon  entering  subroutine  FREEF2,  several  variables  are  initialized.  The  sub¬ 
routine  then  branches  using  an  IF-THEN-ELSEIF  block  that  calculates  the  values  for 
the  arrays  FORCD  and  STRS  depending  on  the  value  of  IVOL(l).  The  procedure  in 
calculating  the  externzJ  forces  for  a  nonreflecting  boundary  is  similar  to  that  used  for  a 
pressure  line.  The  primary  difference  is  that  the  velocity-dependent  viscous  stresses  (or 
pressures)  of  the  nonreflecting  boundary  ate  continually  changing  during  the  solution, 
while  the  pressures  of  the  pressure  load  line  are  dictated  by  the  input  information.  For 
each  successive  pair  of  node  numbers,  or  segments,  along  the  boundary,  the  viscous 
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forces  are  computed  from  the  horizontal  and  vertical  velocities  of  the  two  nodes  and 
the  dilatational  and  shear  impedances  of  the  adjacent  element  within  a  DO  loop.  At 
the  beginning  of  the  DO  loop,  the  horizontal  and  vertical  distances  (SS  and  CC)  be¬ 
tween  the  current  two  nodes  are  calculated.  The  positive  sign  convention  is  shown  in 
Figure  3.3.  The  horizontal  and  vertical  velocities  of  the  two  nodes  are  retrieved  from 
the  array  V.  Because  the  array  BIMP  was  paissed  to  FREEF2  from  the  same  location 
as  KPRES,  the  index  of  BIMP  is  increased  by  NDNOD,  the  number  of  nodes  on  the 
load  line.  Equations  3.32  and  3.33  are  used  to  solve  for  the  viscous  forces.  If  the  non¬ 
reflecting  boundary  is  for  an  axisymmetric  problem,  the  viscous  forces  are  multiplied 
by  the  average  distance  the  two  nodes  are  from  the  axis  of  symmetry.  The  resulting 
viscous  forces  are  finally  stored  in  the  appropriate  locations  in  the  external  force  array 
FORCD. 

The  equations  for  computing  two  of  the  global  variables  for  a  load  line  have  been 
modified -by  the  authors.  The  external  work  done  by  the  viscous  forces  is  incrementally 

added  for  each  segment  along  the  load  line. 

NDNOD 

X)  +  (3-34) 

»=i 


where,  Ax  =  displacement  in  the  z-direction  during  a  a  time  step 

Ay  =  displacement  in  the  y-direction  during  a  a  time  step 
The  equation  for  the  change  in  impulse  on  the  load  line  is 


A/  =  6t, 


NDNOD 


(3.35) 


where, 


6t  =  the  time  step 
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3.4  Revised  Input  Instructions  for  SAMSON2 

The  nonreflecting  boundary  is  treated  as  another  load  line  in  the  input  instruc¬ 
tions.  The  new  updated  load  line  control  data  card  (card  lOA)  is  given  in  Appendix 
A.  Card  lOA  is  the  only  card  in  the  SAMSON2  input  instructions  that  heis  been  mod¬ 
ified.  The  only  modifications  to  the  load  line  card  are  the  extra  two  options  for  the 
load  line  type  (IVOL(l)).  For  a  nonreflecting  boundary  in  a  axisymmetric  problem, 
IVOL(l)  must  equal  three;  and  for  an  plane  problem,  IVOL(l)  must  equal  four.  The 
only  requirement  for  the  use  of  the  nonreflecting  boundary  is  that  the  nodes  on  the 
nonreflecting  boundary  must  be  inputted  in  a  counterclockwise  fashion.  An  example 
is  seen  in  Figure  3.4.  Assume  the  surface  along  the  left  and  bottom  edges  of  the  figure 
are  to  be  nonreflective.  The  nodes  would  be  inputted  as  3,  4,  5,  ...,  80,  90.  Also,  the 
node  generation  option  is  still  applicable  for  the  nonreflecting  boundary  if  the  difference 
between  successive  node  numbers  is  constant. 

Example  input  files  utilizing  the  nonreflecting  boundary  are  located  in  Appendix 
D.  The  nonreflecting  boundary  input  instructions  are  the  last  lines  in  rdl  of  the  input 
files.  The  first  set  of  load  line  instructions  (load  line  1)  is  the  actual  pressure  or 
displacement  function  applied  to  the  system. 

For  4NQ  element  problems  involving  the  AFWL  “engineering"  materird  model, 
Young’s  modulus  is  used  only  for  the  antihourglassing  damping  forces.  The  Poisson’s 
ratio  on  card  5C-1  is  not  used  at  all.  The  values  of  Young’s  modulus  and  Poisson’s  ratio 
should  be  such  that  they  predict  the  approximate  dilatational  and  shear  wave  speeds 
of  the  medium  during  a  solution.  In  the  following  chapter,  the  results  show  that  the 
best  overall  approximations  for  Young’s  modulus  and  Poisson’s  ratio  correspond  to 


Figure  3.4  Order  of  Nodes  Along  a  Nonreflecting  Boundary 


the  values  of  the  bulk  modulus  and  Poisson’s  ratio  for  the  first  loading  segment.  The 
antihourglass  damping  coefficient,  fig,  should  be  used  to  regulate  the  magnitude  of  the 
antihourglass  damping  forces  for  a  given  Young’s  modulus. 
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Chapter  4 

THE  EVALUATION  OF  THE  NONREFLECTING  BOUNDARY 

Once  the  algorithm  for  the  nonreflecting  boundary  was  implemented  into  the 
SAMSON2  code,  a  series  of  tests  was  performed  to  check  the  elTectiveness  of  the  bound¬ 
ary.  Each  individual  test  consisted  of  three  solutions.  The  first  solution  was  obtained 
for  a  large,  or  extended,  mesh.  The  large  mesh  solution  would  simulate  an  infinite 
medium.  A  small  portion  of  the  large  mesh  surrounding  the  loading  area  was  then 
isolated  and  analyzed  under  two  different  sets  of  boundary  conditions.  The  boundary 
conditions  were  defined  along  the  line  that  truncated  the  larger  mesh.  The  second  solu¬ 
tion  was  the  response  of  the  small,  or  truncated,  mesh  with  the  nonreflecting  boundary 
active  along  this  line.  The  comparison  of  these  two  solutions  would  demonstrate  if  the 
nonreflecting  boundary  were  absorbing  the  radiant  wave  energy  and,  thus,  simulating 
an  infinite  medium.  The  third  or  final  solution  was  the  analysis  for  the  same  small 
mesh  but  with  a  fixed/free  boundary  condition.  The  small  mesh  solution  without  the 
nonreflecting  boundary  gave  a  qualitative  measure  of  the  effectiveness  of  the  nonreflect¬ 
ing  boundary  and  showed  what  can  occur  when  the  reflections  of  the  stress  waves  are 
allowed  to  return  and  contaminate  the  analysis  in  the  interior  of  the  mesh.  The  total 
solution  time  for  all  three  solutions  was  limited  to  the  apprcodmate  time  the  first  stress 
wave  reflection  from  the  boundary  of  the  large  mesh  returned  to  the  an^dysis  area  for 
the  three  tests.  This  time  would  determine  when  the  large  mesh  solution  would  cease 
to  behave  as  an  infinite  medium. 
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The  three  solutions  obtained  for  each  test  may  be  inaccurate  due  to  the  probable 
existence  of  logical  errors  in  the  SAMSON2  code.  This  inaccuracy,  however,  is  unlikely 
to  pose  problems  in  the  evaluation  of  the  nonreflecting  boundary.  The  objective  for 
each  test  is  to  compare  the  results  of  three  solutions.  Since  the  nonreflecting  boundary 
depends  solely  upon  parameters  near  the  boundary,  it  does  not  affect  the  behavior  of 
the  interior  of  the  mesh.  The  same  errors  that  may  be  generated  in  one  solution  will  be 
consistently  repeated  in  the  other  two.  Therefore,  the  effectiveness  of  the  nonreflecting 
boundary  will  be  independent  of  the  potential  errors  generated  in  the  three  solutions. 

There  are  four  possible  results  that  can  occur  in  each  test.  The  following  are 
conclusions  which  can  be  obtained  from  each  result: 

1.  If  all  three  of  the  solutions  are  similar; 

the  radiating  stress  waves  in  the  small  mesh  without  the  nonreflecting  boundary 
have  not  had  time  to  propagate  to  the  boundary  and  return  to  the  point  under 
analysis.  Therefore,  the  total  solution  time  should  be  extended  to  allow  sufflcient 
time  for  the  reflections  to  return. 

2.  If  all  three  of  the  solutions  are  not  similar; 

the  nonreflecting  boundary  is  only  partially  effective  in  the  absorption  of  the 
stress  waves.  This  is  an  unsatisfactory  condition. 

3.  If  only  the  small  mesh  solutions  with  and  without  the  nonreflecting  boundary 


are  similar; 
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the  nonreflecting  boundary  is  completely  ineffective.  This  is  the  worst  possible 
result  that  can  occur. 

4.  If  only  the  solution  of  the  small  mesh  with  the  nonreflecting  boundary  and  the 
large  mesh  solution  are  similar, - 

the  nonreflecting  boundary  is  satisfactorily  absorbing  the  radiating  stress  waves. 
The  test  is  successful. 

One  other  result  might  occur  that  is  unrelated  to  the  c^dculations  within  the  SAMSON2 
code.  If  all  three  solutions  do  not  coincide  up  to  a  certain  point,  then  an  error  exists 
in  one  or  more  of  the  input  files.  This  was  a  check  to  verify  that  each  input  file  was 
correct.  If  all  three  of  the  input  files  were  correct,  the  files  are  said  to  be  compatible. 

The  evaluation  of  the  nonreflecting  boundary  is  separated  into  two  categories 
of  problems.  First,  a  simple  one-dimensional  wave  propagation  problem  is  thoroughly 
investigated.  Many  facets  of  the  SAMSON2  code,  including  both  materisd  and  element 
properties,  were  varied  for  a  bar  subjected  to  a  dilatational  wave  and  a  shear  wave. 
This  study  provided  simple  tests  for  the  effectiveness  of  the  nonreflecting  boundary 
for  several  parameters.  Also,  the  various  tests  in  the  one-dimensional  wave  propaga¬ 
tion  study  was  used  to  check  the  formulation  of  the  nonreflecting  boundary  within  the 
SAMSON2  code.  Once  the  nonreflecting  boundary  appeared  to  be  functioning  cor¬ 
rectly,  a  much  more  realistic  and  complex  problem  was  analyzed.  The  second  part  of 
the  evaluation  involved  a  two-dimensional  axisymmetric  problem  that  provided  tests 
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to  check  the  effectiveness  of  the  nonreflecting  boundary  subject  to  stress  waves  with 
varying  angles  of  incidence  and  material  properties. 

The  original  scope  of  the  research  was  to  include  the  solution  of  a  typical  soil- 
structure  interaction  (SSI)  problem  that  has  been  analyzed  by  the  personnel  at  the 
AFWL.  Unfortunately,  the  slideline  interface  routine  within  the  authors’  version  of 
SAMSON2  did  not  function  properly.  The  time  to  correct  the  errors  in  the  slideline 
routine  would  not  justify  the  need  to  include  a  slideline  interface  in  the  test  prob¬ 
lems.  The  nonreflecting  boundary  is  independent  of  the  interior  properties  of  the 
mesh.  Therefore,  the  slideline  interface  would  not  affect  the  effectiveness  of  the  nonre¬ 
flecting  boundary  and,  thus,  its  inclusion  in  the  two-dimensional  problems  would  not 
be  essential  in  this  study. 

4.1  One- Dimensional  Wave  Propagation 

The  initial  tests  involving  the  nonreflecting  boundary  consisted  of  examining 
various  sample  problems  that  were  modiflcations  of  the  one-dimensional  wave  prop¬ 
agation  problem  in  Appendix  A  of  the  SAMSON2  user’s  manual  [12].  The  purposes 
of  utilizing  this  particular  problem  are  twofold.  The  computer  solutions  that  exist 
in  the  user’s  manual  provided  the  opportunity  to  check  that  the  present  version  of 
SAMSON2,  transferred  to  the  Washington  State  University’s  IBM  3090,  performs  cor¬ 
rectly.  Second,  a  comparison  of  the  original  SAMSON2  solutions  in  the  users  manual 
to  the  modified  version  of  SAMSON2  solutions  affirmed  that  the  implementation  of 
the  nonreflecting  boundary  into  the  code  did  not  affect  the  correct  results. 
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Two  bars  were  used  exclusively  throughout  the  one-dimensional  study.  A  long 
bar,  which  will  be  referred  to  as  a  large  mesh,  was  40-cm  long  by  4-cm  high  with  a 
thickness  of  1-cm.  The  right  end  of  the  long  bar  was  free  to  translate.  A  short  bar, 
accordingly  referred  to  as  a  small  mesh,  had  the  same  cross-sectional  dimensions  as 
the  long  bar  while  the  length  was  shortened  to  6-cm.  In  one  case,  the  nonreflecting 
boundary  was  applied  to  the  right  end  of  the  short  bar;  in  the  other  ceise,  the  right  end 
was  free.  Each  bar  is  perfectly  elastic  and  has  no  internal  damping.  Poisson’s  ratio  is 
zero  except  when  noted  and  Young’s  modulus  and  the  mass  density  are  each  assigned 
the  value  of  one. 

For  all  the  wave  propagation  tests,  a  displacement  history  function  described  in 
Figure  4.1  wxs  applied  to  the  left  end  of  each  of  the  bats.  When  a  bar  is  subjected  to 
a  dilatational  wave,  all  nodes  ate  free  to  translate  in  both  the  horitontal  and  vertical 
direction  relative  to  the  longest  axis  of  the  bar.  The  displacement  function  is  then 
applied  normal  to  the  nodes  on  the  left  face  of  the  bar.  When  a  bar  is  subjected  to  a 
shear  wave,  all  nodes  are  free  to  translate  in  the  vertical  direction  but  are  fixed  in  the 
horizontal  direction.  This  restricts  the  bar  to  deform  only  in  shear  and  greatly  reduces 
the  flexural  deformation.  The  displacement  function  is  then  applied  tangentirdly  to  the 
nodes  on  the  left  face  of  the  bar. 

All  of  the  various  continuum  elements  in  the  SAMSON2  library  were  tested. 
However,  only  the  four-node  quadrilateral  (4NQ)  and  the  eight-node  quadrilateral 
(8NQ)  element  results  are  presented.  The  4NQ  and  8NQ  are  the  principal  elements 
used  by  the  personnel  at  the  AFWL.  The  remaining  triangular  elements  are  mostly  used 
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Figure  4.1  Displacement  Function 

in  a  transition  between  the  quadrilateral  elements  and  were  not  considered  to  be  im¬ 
portant  in  this  study.  The  results  of  the  one-dimensional  problem  using  the  triangular 
elements  were  similar  -  but  slightly  less  accurate  -  to  the  results  using  the  quadrilateral 
elements.  The  SAMSON2  input  files  for  a  selected  number  of  wave  propagation  tests 
are  listed  in  Appendix  D. 

The  tests  that  involved  the  4NQ  element  used  the  two  discretized  meshes  shown 
in  Hgute  4.2.  The  large  mesh  consists  of  63  nodes  smd  40  elements.  The  small  mesh 
has  12  nodes  and  6  elements.  Each  element  is  2-cm  by  2-cm  with  a  thickness  of  1- 
cm.  The  horizontal  and  vertical  displacement  of  a  node  located  2-cm  from  the  left  and 
2-cm  from  the  bottom  of  the  bar  (x  =  2-cm,  y  =  2-  cm)  are  compared  in  the  three 


solutions.  The  normal  and  shear  stresses  are  also  evaluated  at  the  center  of  an  element 


40-cm 


42 


LARGE  MESH 
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Direction  of  Wave  Propagation 
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Figure  4.2  4NQ  Discretization  of  One-Dimensional  Bars 
at  a  point  3-cm  from  the  left  and  l-tm  from  the  bottom  of  the  bar  (x  —  3-cm,  y  = 
1-cm).  The  positions  of  both  points  in  the  two  discretized  bars  are  marked  as  dots  in 
the  figure. 

The  first  test  was  a  bar  consisting  of  4NQ  plane  continuum  elements  subjected 
to  a  dilatational  stress  wave.  The  graphical  results  for  the  displacement  and  the  normal 
stress  versus  time  are  shown  in  Figures  4.3  and  4.4.  The  large  mesh  solution,  simulating 
an  infinitely  long  bar,  is  plotted  as  a  solid  line.  The  short  dashed  line  represents  the 
small  mesh  solution  with  the  nonreflecting  boundary  applied  at  the  right  end  of  the  bar. 
If  the  nonreflecting  boundary  is  working  properly,  the  solid  line  and  the  short-dashed 
line  should  be  superimposed.  Finally,  the  long-dashed  line  represents  the  small  mesh 
with  a  free  right  end  and  demonstrates  the  effective  absorption  of  the  nonreflecting 
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boundary.  Note  that  all  three  curves  coincide  up  to  a  time  of  about  8-sec  indicating 
that  the  input  files  are  compatible.  The  two  figures  show  that  there  is  almost  complete 
absorption  of  the  dilatational  stress  wave. 

The  next  test  was  similar  to  the  previous  test  except  that  the  bars  are  now 
subjected  to  a  shear  wave  instead  of  a  dilatational  wave.  The  graphical  results  are 
presented  in  Figures  4.5  and  4.6.  Once  ag2un,  there  is  nearly  perfect  absorption  of  the 
shear  wave.  The  slight  discrepancy  in  the  shear  wave  solutions  can  be  accounted  for 
by  small  amounts  of  flexural  deformations  that  occur  in  the  bars. 

At  this  point  in  the  testing,  several  conclusions  could  be  made.  The  edgorithm 
of  the  nonreflecting  boundary  implemented  into  the  subroutine  FREEF2  is  correct 
for  cill  plane  continuum  elements.  The  procedure  for  finding  the  elements  along  the 
nonreflecting  boundary  for  the  4NQ  continuum  element  and  the  equations  for  the 
dilatational  and  shear  impedences  written  into  the  subroutine  READPV  are  correct. 
There  were  no  errors  in  the  passage  of  the  various  arrays  added  to  the  argument  lists 
of  the  modified  subroutines  which  suggest  tha.t  the  storage  pointers  for  the  arrays  have 
been  correctly  indexed. 

To  insure  the  complete  accuracy  for  the  algorithm  of  the  nonreflecting  boundary 
in  the  subroutine  FREEF2,  tests  were  conducted  with  the  4NQ  axisymmetric  contin¬ 
uum  element.  Bars  parallel  and  perpendicular  to  the  axis  of  symmetry  were  again 
subjected  to  a  dilatational  and  a  shear  wave.  For  the  bars  parallel  to  the  ws  of 
symmmetry,  identical  solutions  to  those  calculated  in  the  previous  plane  case  were 


observed. 
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Figure  4.3  Horiiontal  Displacement  vs.  Time  at  x  =  2-cm, 
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Figure  4.5  Vertical  Displacement  vs.  Time  at  x 
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For  the  bars  perpendicular  to  the  axis  of  symmetry,  the  geometry  of  the  ax- 
isymmetric  problem  had  to  be  altered.  White,  Valliappan,  and  Lee  [20]  stated  that 
the  efficiency  of  the  viscous  nonreflecting  boundary  was  dependent  upon  the  ratio  of 
the  distance  the  boundary  w^«  from  the  axis  of  symmetry  to  the  wavelength  of  the 
incident  stress  wave.  Recall  that  for  an  R/A  ratio  greater  than  0.5,  the  nonreflecting 
boundary  should  be  effective.  The  sinusoidal  displacement  function  applied  to  the  left 
end  of  the  bar  has  an  angular  frequency  of  0.1963  rad/sec  (f  =  0.0312  Hz).  The  wave 
speed  of  the  dilatational  wave  through  the  bar  is  1.0  cm/sec.  It  can  be  easily  proven 
that  the  wavelength  of  the  dilatational  wave  is  32  cm.  Two  tests  were  conducted  using 
R/A  values  of  0.5  and  2.0,  respectively.  For  an  R/A  vaJue  of  0.5,  the  two  bars  were 
placed  such  that  the  left  ends  of  the  bars  were  10-cm  from  the  axis  of  symmetry.  The 
right  end  of  the  short  bar  would  then  be  16-cm  from  the  axis  of  symmetry.  Similarly 
for  an  R/A  value  of  2.0,  the  left  ends  of  the  bars  are  located  at  x  =  58-cm.  Only  the 
dilatational  wave  graphical  results,  shown  in  Figures  4.7  •  4.10,  are  given  for  the  two 
values  of  R/A.  Consistent  with  reference  [20],  u  the  ratio  of  R/A  is  increased,  the 
efficiency  of  the  nonreflecting  boundary  does,  in  fact,  increase.  Note  that  for  R/A  = 
0.5,  the  results  do  not  correspond  as  well  as  those  for  R/A  =  2.0.  Therefore,  to  insure 
adequate  results,  the  limiting  R/A  value  of  0.5,  should  be  increased  to  2.0. 

With  the  successful  completion  of  the  axisymmetric  problem,  the  formulation  of 
the  nonreflecting  boundary  within  the  subroutine  FREEF2  was  determined  to  be  error- 
free.  The  only  remaining  untested  portion  of  the  code  is  the  procedure  in  subroutine 
READPV  that  retrieves  the  element  numbers  along  the  nonreflecting  boundary.  The 
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Figure  4.7  Horizontal  Displacement  vs.  Time  at  x  =  12-cm,  y  ==  2-cm  for 
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Figure  4.8  a,,  vs.  Time  at  x  =  13-cm,  y  =  1-cin  for  R/A  =  0.5 
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Figure  4.9  Horizontal  Displacement  vs.  Time  at  x  =  60-cm,  y  =  2  cm  for  R/A 
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Figure  4.11  8NQ  Discretization  of  One- Dimensional  Bars 
triangular  elements  and  the  eight-node  quadrilateral  elements  still  required  verification. 
For  each  of  these  elements,  only  a  dilatational  stress  wave  test  was  used. 

The  40-cm  bar  and  6-cm  bar  were  now  discretized  using  2-cm  by  4-cm  8NQ  plane 
isoparametric  elements.  The  meshes  are  shown  in  Figure  4.11.  The  long  bar  consists 
of  20  elements  and  103  nodes  while  the  short  bar  has  3  elements  and  18  nodes.  All  of 
the  material  properties  remain  unchanged.  The  integration  order  of  the  8NQ  elements 
was  two.  Figures  4.12  and  4.13  give  the  graphical  results  of  the  dilatational  stress  wave 
tests.  The  displacement  versus  time  graph  for  the  plane  8NQ  is  very  similar  to  that  of 
the  plane  4NQ  element.  The  graph  showing  the  normal  stress  versus  time,  although 
exhibiting  very  good  results,  is  more  erratic  in  comparison  to  the  corresponding  graph 
of  the  plane  4NQ  element. 
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Figure  4.12  Hoiitontal  Displacement  vs.  Time  at  x  =  2-cm,  y  =  2-cm 
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Figure  4.13  a,,  vs.  Time  at  x  =  2.423  cm,  y  =  0.423  cm 
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Similar  tests  were  conducted  for  the  three-node,  five-node  and  six-node  plane 
triangular  elements.  Although  the  results  have  not  been  included  in  this  report,  all 
three  triangular  elements  displayed  the  same  behavior  as  the  4NQ  and  8NQ  elements 
except  for  slightly  less  accuracy.  At  the  conclusion  of  these  tests,  all  of  the  modifications 
to  the  SAMSON2  code  were  believed  to  be  correct. 

The  purpose  of  the  remaining  tests  in  this  research  was  to  assess  the  versatility 
of  the  nonreflecting  boundary  under  varying  materi^d  properties.  Poisson’s  ratio,  i/,  of 
the  medium  is  known  to  affect  the  effectiveness  of  the  nonreflecting  boundary.  Two 
tests  were  conducted  employing  the  plane  4NQ  discretized  bars  for  two  different  values 
of  Poisson’s  ratio:  v  =  0.20  and  v  =  0.45,  A  value  of  0.45  would  be  an  extreme  value 
of  Poisson’s  ratio  (i/  =  0.5  would  imply  an  incompressible  material).  The  bars  were 
again  subjected  to  a  dilatationai  wave  and  a  shear  wave,  respectively.  Figures  4.14 
-  4.17  show  the  solutions  associated  with  v  =  0.2.  There  is  a  small  deviation  in  the 
graphs  shown  in  the  four  figures  between  the  large  mesh  solution  and  the  nonreflecting 
boundary  solution.  This  deviation  becomes  greater  for  v  =  0.45  as  is  evident  by  the 
graphs  shown  in  Figures  4.18  -  4.21.  The  results  from  these  tests  support  the  evidence 
of  Cohen  and  Jennings  [3]  that  the  efficiency  of  the  nonreflecting  boundary  decreases 
as  Poisson’s  ratio  increases.  Furthermore,  the  degree  of  error  demonstrated  in  the 
two  tests  is  appreciable  only  when  Poisson’s  ratio  tends  to  its  maximum  value  of  0.5. 
For  Poisson’s  ratio  in  the  range  of  0  to  0.4,  the  nonreflecting  boundary  performed 
satisfactorily. 

The  final  test  of  the  nonreflecting  boundary  involving  one-dimensionaJ  wave 
propagation  used  the  plane  4NQ  discretized  bar  consisting  of  two  different  materials 


Figure  4.14  Horizontal  Displacement  vs.  Time  at  x  =  2-cm,  y  =  2-cin  fo 
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Figure  4.16  Vertical  Displacement  vs.  Time  at  x  =  2-cm,  y  —  2-cm  for  i/ 
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Figure  4.18  Horizontal  Displacement  vs.  Time  at  x  =  2-cm, 
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Figure  4.20  Vertical  Displacement  vs.  Time  at  x  =  2-cm, 
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Figure  4.21  a,,,  vs.  Time  at  x  =  3-cm,  y  =  1-cm  for  i/  =  0.45 
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(i.e.,  an  anisotropic  medium).  The  lower  half  of  both  bars  had  the  same  material 
properties  as  before  with  u  =  0.0.  Only  Young’s  modulus  for  the  upper  halves  of  the 
bars  differed  from  the  lower  halves  of  the  bars.  Young’s  modulus  for  the  upper  halves  of 
the  bars  was  4.0  dyne/cm^.  Therefore,  the  dilatational  wave  speed  and  the  shear  wave 
speed  for  the  upper  half  is  double  that  for  the  lower  half.  The  bars  were  again  subjected 
to  a  dilatational  wave  and  a  shear  wave.  Figures  4.22  -  4.25  show  the  graphical  results 
for  these  tests.  The  nonreflecting  boundary  absorbed  both  stress  waves  very  effectively. 
Although,  this  test  may  be  oversimplified  to  conclude  that  the  nonreflecting  boundary 
can  adequately  absorb  stress  waves  in  a  two-dimensional  space,  it  does,  however,  give 
an  indication  of  its  capability. 

4.2  Two-Dimensional  Wave  Propagation 

The  second  category  of  the  evaluation  of  the  nonreflecting  boundary  entails 
the  use  of  a  two-dimensional  half-space.  Within  a  two-dimensional  half-space  the 
dilatational  waves  and  the  shear  waves  can  impinge  simultaneously  on  the  nonreflecting 
boundary  at  different  angles  of  incidence.  The  angles  of  incidence  of  the  stress  waves  are 
unknown  during  the  analysis  of  the  system.  The  stress  and  displacement  distributions 
in  a  mesh  are  results  of  the  combined  effects  of  both  the  dilatational  wave  and  the  shear 
wave.  It  is  difficult  to  distinguish  the  contribution  of  an  individual  stress  wave.  The 
response  of  the  system  due  to  these  radiant  stress  waves  can  become  quite  complicated. 
In  the  evaluation  of  the  nonrefiecting  boundary,  therefore,  the  emphasis  is  not  on  the 
actual  response  of  the  system  (unless,  of  course,  the  results  are  obviously  incorrect) 
but  on  the  comparison  of  the  three  solutions  during  each  test. 
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Figure  4.22  Horizontal  Displacement  vs.  Time  at  x  =  2-cm,  y  -  2-cm 
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Figure  4.24  Vertical  Displacement  vs.  Time  at  x  =  2-cm,  y  =  2-cm 
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The  large  mesh  that  will  be  used  in  the  two-dimensional  study  is  shown  in 
Figure  4.26.  It  is  composed  of  6-in  by  6-in  4NQ  axisymmetric  continuum  elements. 
The  dimensions  of  the  large  mesh  are  90-in  by  90-in.  The  left  edge  and  the  right 
edge  are  fixed  in  the  horizontal  direction  and  free  in  the  vertical  direction.  For  aii 
axisymmetric  case,  the  left  edge  will  lie  along  the  axis  of  symmetry.  The  base  is  fixed 
in  the  vertical  direction  and  free  in  the  horizontal.  The  fixity  conditions  of  the  three 
edges  result  in  the  bottom  two  corner  nodes  being  fixed  in  both  directions.  The  loading 
will  be  applied  to  the  top  edge  of  the  mesh.  The  dotted  line  in  the  upper  left  section 
of  the  mesh  marks  the  truncation  line  that  defines  the  boundary  of  the  corresponding 
small  mesh.  The  small  mesh  is  shown  in  Figure  4.27.  For  one  solution,  the  nonreflecting 
boundary  is  applied  along  this  truncation  line.  For  the  next  solution,  the  small  rnesh 
has  the  exact  fixity  conditions  as  that  of  the  large  mesh. 

Each  two-dimensional  wave  propagation  test  involves  three  field  variables  versus 
time  at  a  certain  point:  the  displacement  of  a  node  located  12-in  from  the  left  edge 
and  24-in  from  the  top  edge  and  the  normal  stress,  o-yy,  and  the  shear  stress,  (Txy,  at 
a  point  15-in  from  the  left  edge  and  21-in  from  the  top  edge.  The  locations  of  these 
two  points  are  shown  in  Figure  4.27.  However,  the  results  from  these  tests  can  be 
misleading.  The  analysis  of  these  field  variables  are  indicative  of  the  response  of  the 
mesh  only  in  the  vicinity  of  these  points  of  observation.  While  the  graphical  results 
may  show  that  the  nonreflecting  boundary  is  effective,  areas  further  away  from  these 
observed  points  may  show  that  the  nonreflecting  boundary  is  completely  ineffective. 
Instead  of  checking  each  point  in  the  mesh  over  time  to  prove  conclusively  that  the 


nonreflecting  boundary  was  performing  well,  the  vertical  displacement  of  each  node  in 
the  mesh  was  examined  at  a  certain  time  for  each  solution. 

Of  the  three  field  variables,  the  displacement  plots  are  relatively  unimportant 
in  comparison  to  the  normal  stress  and  shear  stress  plots.  There  will  be  a  “residual'’ 
displacement  at  the  completion  of  a  problem  utilizing  a  nonreflecting  boundary.  The 
displacement  of  the  nodes  in  the  two  meshes  are  given  simply  to  provide  a  general 
indication  of  the  effectiveness  of  the  nonreflecting  boundary. 

The  limiting  solution  time  for  the  three  field  variables  was  the  approximate 
time  the  reflections  returned  from  the  boundary  of  the  large  mesh.  Unfortunately,  it 
was  impossible  to  tell  when  these  reflections  actually  did  return.  The  distance  from 
the  point  of  excitation  to  the  boundary  of  the  large  mesh  back  to  the  location  of  the 
field  variables  w«  roughly  three  times  the  distance  for  the  same  route  in  the  small 
mesh.  Assuming  that  the  stress  waves  travelled  at  the  same  speed  in  both  meshes,  the 
time  for  the  stress  wave  to  return  from  the  boundary  of  the  large  mesh  would  be  three 
times  the  time  for  the  stress  wave  to  return  from  the  boundary  of  the  small  mesh. 
Therefore,  from  the  time  the  three  solutions  began  to  diverge  (the  time  the  reflections 
have  returned  from  the  boundary  of  the  small  mesh),  the  limiting  solution  time  was 
computed  to  be  three  times  the  time  at  this  divergent  point. 

The  material  law  used  exclusively  for  the  two-dimensional  problems  was  the 
AFWL  “engineering”  model.  Two  different  materials  were  modelled:  Yuma  soil  and 
fiber  reinforced  concrete.  The  mateti^d  parameters  of  the  Yuma  soil  are  listed  in  Figure 
4.28  along  with  the  hydrostatic  pressure  versus  volumetric  strain  curve  and  the  yield 
surface  curve.  The  same  items  are  shown  in  Figure  4.29  for  fiber  reinforced  concrete. 


Hydrostatic  Pressure  {L US/ 1 ) 


Mass  Density  =  0.173E'3  lbs  ■  sec' /in^  Poisson’s  Ratio  =  0.38 

Young’s  Modulus  =  0.1904E5  Ibsfsec^  Fraction  of  Critical  Damping  =  0. 

Figure  4.28  Material  Properties  of  Yuma  Soil 


yjj  {LUSjlN^)  Hydrostatic  Pressure  (Liy5//yV*) 


Mass  Density  =  0.2516E-3  lbs  •  sec^ /in*  Poisson’s  Ratio  =  0.24 

Young’s  Modulus  =  0.1755E7  Ibs/sec^  Fraction  of  Critical  Damping  =  0 

Figure  4.29  Material  Properties  of  Fiber-Reinforced  Concrete 


Note  that  the  values  of  Young’s  modulus  and  Poisson’s  ratio,  which  determine  the 
dilatational  and  shear  impedences  of  the  material,  correspond  to  the  bulk  modulus 
and  Poisson’s  ratio  of  the  first  loading  segment  for  both  materials. 

For  all  of  the  two-dimensional  wave  propagation  problems,  the  large  mesh  and 
the  small  mesh  were  isotropic.  It  Wcis  impractical  to  analyze  an  anisotropic  medium. 
If  the  large  mesh  were  anisotropic,  the  radiating  stress  waves  would  both  refract  and 
reflect  at  the  interface  for  a  change  in  material  properties.  The  stress  wave  reflections 
from  this  interface  may  return  to  the  area  under  analysis  prematurely.  The  large  mesh 
would  then  cease  to  behave  as  an  infinite  medium.  Therefore,  it  would  not  be  applicable 
to  compare  the  nonreflecting  boundary  solution  to  the  large  mesh  solution,  and  the 
test  would  be  invalid. 

Several  tests  were  initially  conducted  employing  Yuma  soil.  The  performance 
of  the  nonreflecting  boundary  was  analyzed  under  a  variety  of  test  conditions.  Unless 
specified  otherwise,  the  loading  function  used  will  be  a  Brode  nuclear  airblast  curve 
with  a  peak  pressure  of  50  ksi.  The  airblast  curve  is  shown  in  Figure  4.30.  Also,  the 
two  meshes  will  consist  of  4NQ  axisymmetric  elements  for  which  the  axis  of  symmetry 
will  lie  along  the  left  edge  of  the  mesh. 

For  both  Yuma  soil  and  concrete,  the  parameters  of  the  AFWL  “engineering” 
material  law  describe  a  very  nonlinear  material.  Together  with  the  large  magnitude 
loading  condition,  the  test  cases  examined  in  this  two-dimensional  study  are  extreme 
tests  of  the  nonreflecting  boundary.  The  effective  limit  of  the  nonreflecting  boundary 
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Figure  4.30  Erode  Nuclear  Airblast  Curve 
was  sought.  For  problems  that  involved  an  elastic  or  near  elastic  material,  the  nonre¬ 
flecting  boundary  performed  very  well  regardless  of  the  type  of  material  law.  However, 
the  results  have  been  omitted  in  this  report. 

The  purpose  of  the  first  two  tests  involving  two-dimensional  wave  propagation 
was  to  examine  the  effect  of  loading  applied  tangentially  along  the  nonreflecting  bound¬ 
ary  line.  For  the  first  test,  the  airblast  load  was  applied  vertically  downward  along  the 
top  line  of  the  small  mesh  from  the  euds  of  symmetry  to  the  top  right  node  of  the 
small  mesh.  Thus,  the  right  end  of  this  load  line  is  in  contact  with  the  nonreflecting 
boundary.  Because  the  load  line  length  for  both  meshes  must  be  equal,  the  airblast 
load  was  applied  from  the  axis  of  symmetry  of  the  large  mesh  to  the  node  located  36-in 
to  the  right. 


The  graphical  results  are  presented  in  Figures  4.31  -  4.33.  Note  that  the  ti.Tie 
of  divergence  of  the  three  solutions  is  approximately  0.0015  sec.  The  termination 
time  for  the  solutions  is  0.004  sec.  Because  the  point  of  divergence  is  greater  than 
1/3  the  termination  time,  the  large  mesh  simulates  an  infinite  medium  during  the 
entire  solution  period.  The  displacement  versus  time  and  normal  stress  versus  time 
graphs  show  a  good  correlation  between  the  large  mesh  solution  and  the  nonreflecting 
boundary  solution.  However,  the  plot  of  the  shear  stress  versus  time  exhibits  more 
discrepancy  than  the  previous  two  graphs. 

The  second  test  was  identical  to  the  previous  test  except  that  the  length  of  the 
airbiast  load  line  was  shortened.  For  both  meshes,  the  load  line  now  extends  from  the 
axis  of  symmetry  for  a  distance  of  18-in.  The  same  three  parameters  versus  time,  as 
were  analyzed  in  the  first  test,  are  examined.  The  graphical  results  of  the  test  are 
shown  in  Figures  4.34  -  4.36.  In  adl  three  cases,  the  nonreflecting  boundary  solution 
was  very  comparable  with  the  large  mesh  solution. 

In  order  to  investigate  the  response  of  the  entire  mesh,  the  vertical  displacement 
distribution  of  all  nodes  were  plottea  at  time  =  0.004  sec.  Figure  4.37  shows  the  vertical 
displacement  of  the  nodes  in  the  upper  left  section  of  the  large  mesh.  Likewise,  Figures 
4.38  and  4.39  show  the  vertical  displacement  of  all  the  nodes  in  the  small  mesh.  To 
clarify  the  orientation  of  the  plots  in  the  figures,  the  bottom  left  corner  of  the  plot 
corresponds  to  the  bottom  left  corner  of  the  smaJl  mesh.  The  vertical  displacements 
are  projected  upward  from  each  node  in  the  mesh.  The  large  mesh  displacement 
distribution  is  bounded  by  28.77-in.  while  the  nonreflecting  boundary  displacement 
distribution  is  bounded  by  28.10-in.  The  surfaces  described  by  the  displacements  in 
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Figure  4.31  Vertical  Displacement  vs.  Time.  Extended  Load  Line 
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Figure  4.33  Shear  Stress,  vs.  Time.  Extended  Load  Line 
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Figure  4.34  Vertical  Displacement  vs.  Time  Shortened  Load  Line 


Figure  4.35  Normal  Stress,  Oyy  vs.  Time  Shortened  Load  Line 
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Figure  4.36  Shear  Si..is,  a^u  vs.  Time  Shorteiieci  Load  Line 
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Figure  4.37  Vertical  Dispiaceiueut  Distribution  for  llu-  Large  Mesli 
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Figure  4.40  Impulse  Load 

Figures  4.37  and  4.38  are  very  similar.  The  vertical  displacement  distribution  of  the 
small  mesh  without  the  nonreflecting  boundary  is  bounded  by  17.5d-in  and  exhibits 
a  surface  that  differs  slightly  from  the  previous  two  plots.  The  vertical  displacement 
along  the  bottom  axis  in  Figure  4.39  is  zero  indicating  that  this  axis  is  indeed  fixed  in 
the  vertical  direction. 

The  two  tests  demonstrate  a  limitation  of  the  nonreflecting  boundary.  The  shear 
absorption  of  the  nonreflecting  boundary  tangent  to  the  vertical  loading  in  the  small 


mesh  does  not  adequately  simulate  the  shear  stiffness  of  the  portion  of  the  large  mesh 
that  is  to  the  immediate  right  of  the  vertical  load.  The  remaining  tests  in  this  study, 
will  apply  the  load  along  the  18-in  line  measured  right  from  the  axis  of  symmetry. 
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The  nonreflecting  boundary  ha^  performed  well  for  the  tests  utilizing  a  Erode 
nuclear  airblaist  curve.  The  loading  was  then  changed  to  a  simple  triangular  impulse 
load  shown  in  Figure  4.40  with  a  peak  pressure  of  50  ksi.  All  other  parameters  for 
both  meshes  have  been  kept  the  same.  The  graphical  results  of  the  three  field  variables 
in  Figures  4.41  -  4.43  demonstrate  that  the  form  of  the  loading  does  not  affect  the 
effectiveness  of  the  nonreflecting  boundary.  In  this  case,  the  results  are  better  than  the 
previous  tests  using  the  airbl2ist  curve.  The  large  mesh  solution  and  the  nonreflecting 
boundary  solution  in  the  normal  stress  and  shear  stress  versus  time  graphs  have  similar 
peak  values  of  stress  but  are  slightly  out  of  phase  with  one  another.  The  cause  of  this 
phenomena  is  unknown  to  the  authors  and  several  unsuccessful  attempts  were  made 
to  remedy  the  situation.  The  vertical  displacement  distributions  of  the  three  meshes 
at  time  =  0.004  sec  are  shown  in  Figures  4.44  -  4.46.  The  surfaces  described  by  the 
vertical  displacements  of  the  nodes  in  the  three  meshes  are  similar  to  those  in  Figures 
4.37  -  4.39  but  are  bounded  by  different  displacements. 

The  rem^dning  tests  of  the  nonreflecting  boundary  incorporated  fiber  reinforced 
concrete  instead  of  Yuma  soil.  The  loading,  once  agetin,  is  the  eurblast  curve  with  a 
peak  pressure  of  50  ksi.  The  graphical  results  are  presented  in  Figures  4.47  -  4.49. 
Note  that  the  termination  time  for  the  solutions  is  0.0015  sec  instead  of  0.004  sec  for 
Yuma  soil.  The  stress  wave  speed  through  concrete  is  approximately  three  times  feister 
than  the  stress  wave  speed  through  Yuma  soil.  The  time  of  divergence  for  the  three 
solutions  is  still  greater  than  1/3  the  termination  time.  Therefore,  the  test  is  valid. 
The  nonreflecting  boundary  solution  correlates  well  with  the  large  mesh  solution  for  all 
three  field  variables,  particularly  the  normal  stress.  The  normal  stress  and  the  shear 
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Figure  4.41  Vertical  Displacement  vs.  Time  Impulse  Load 
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Figure  4.42  Normal  Stress,  Oyy  vs.  Time  Impulse  Load 
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Figure  4.43  Shear  Stress,  vs.  Time  Impulse  Load 
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Figure  4.44  Vertical  Displacement  Distribution  for  the  Large  Mesh 
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Figure  4.45  Vertical  Displacement  Distribution  for  the  Small  Mesh  with  a  Nonreflecting  Boundary 
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Figure  4.46  Vertical  Displacement  Distribution  for  the  Small  Mesh  without  a  Nonrenecting  Boundary 
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Figure  4.47  Vertical  Displacement  vs.  Time.  Fiher-lteiiifi 
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Figure  4.50  Vertical  Displacement  Distribution  fi>r  tiie  barge  Mesh 
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Figure  4.51  Vertical  Displacement  Distribution  for  the  Small  Mesh  with  a  N<»/irelle<  ling  Honiidary 
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stress  versus  time  graphs  are  not  as  erratic  as  the  corresponding  graphs  involving  Yuma 
soil.  The  vertical  displacement  distribution  of  the  upper  left  section  of  the  large  mesh 
is  shown  in  Figure  4.50.  The  vertical  displacement  distribution  of  the  small  mesh  with 
the  nonreflectmg  boundary  is  shown  in  Figure  4.51.  The  surfaces  of  these  two  plots 
are  reasonably  close  to  one  another.  This  discrepancy  is  negligible  upon  observation 
of  the  vertical  displacement  distribution  of  the  small  mesh  without  the  nonreflecting 
boundary  illustrated  in  Figure  4.52. 

.\11  of  the  previous  tests  used  loads  with  a  peak  pressure  of  50  ksi.  The  magni¬ 
tude  of  this  peak  pressure  caused  large  inelastic  deformations  to  occur  in  the  meshes. 
The  next  test  employed  the  same  form  of  the  Erode  Airblast  curve  but  with  a  lower 
peak  pressure  of  50  psi.  For  fiber-reinforced  concrete,  the  first  loading  segment  corre¬ 
sponds  to  the  unloading  segment  in  the  hydrostatic  pressure  versus  volumetric  strain 
curve.  At  low  levels  of  stress,  the  concrete  mesh  will  remain  elwtic.  The  first  loading 
segment  for  Yuma  soil,  however,  is  much  different  than  the  unloading  segment.  At 
low  levels  of  stress,  the  soil  will  still  exhibit  inelastic  behavior.  The  graphical  results 
of  the  solutions  involving  Yuma  soil  at  small  magnitude  of  loading  were  extremely 
noisy.  It  was  difficult  to  distinguish  between  the  three  solutions.  Therefore,  only  the 
test  involving  concrete  is  included  in  this  report  for  low  peak  pressure  loadings.  The 
displacement,  normal  stress,  and  shear  stress  versus  time  for  the  selected  field  variables 
were  analyzed.  The  graphical  results  are  presented  in  Figures  4.53  -  4.55.  Once  again, 
the  large  mesh  solutions  were  very  similar  to  the  nonreflecting  boundary  solutions.  The 
displacement  versus  time  graph  shows  more  variance  between  the  three  solutions  while 
a  less  significant  difference  can  be  seen  in  the  normal  stress  and  shear  stress  versus 
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time  graphs.  The  vertical  displacement  distributions  of  the  meshes  were  not  plotted 
due  to  the  similarity  to  the  previous  test.  This  test  is  indicative  of  the  ability  of  the 
nonreflecting  boundary  to  absorb  stress  waves  of  varying  strengths. 
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Figure  4.53  Vertical  Displacement  vs.  Time  Lower  Peak  Pressure 
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Figure  4.S5  Shear  Stress,  vs.  Time  Lower  Peak  Pressure 
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Chapter  5 

CONCLUSIONS  AND  RECOMMENDATIONS 

The  nonreflecting  boundary  that  has  been  implemented  into  the  SAMSON2 
code  has  proven  to  be  very  efTective  in  the  absorption  of  dilatational  stress  waves  and 
shear  stress  waves  under  a  wide  range  of  applications.  There  are,  however,  certain 
limitations  on  the  use  of  the  nonreflecting  boundary.  Some  of  these  limitations  have 
not  been  fully  tested  and  should  requite  further  analyses.  There  are  many  conclusions 
that  can  be  obtained  from  the  results  of  the  several  tests  conducted  in  this  study.  Rec¬ 
ommendations,  in  addition  to  conclusions,  are  presented  in  order  to  provide  instruction 
on  the  proper  and  most  efficient  use  of  the  nonreflecting  boundary. 

The  objective  of  a  nonreflecting  boundary  is  to  minimize  the  computational 
and  storage  costs  for  the  computer  and  the  time  required  by  the  user  to  prepare  an 
input  file.  Because  the  boundary  is  treated  as  a  load  line,  the  additional  computation 
time  and  storage  for  the  input  instructions  for  the  boundary  are  about  the  same  as 
if  another  load  line  were  applied  to  the  system.  Also,  the  input  instructions  for  the 
boundary  are  easily  understood  by  users  familiar  with  the  SAMSON2  code  and  the 
format  of  the  load  line  data  card.  The  introduction  of  the  nonreflecting  boundary  into 
the  SAMSON2  code  required  relatively  few  executable  statements.  The  additional 
execution  time  for  the  nonreflecting  boundary  is  a  result  of  the  modified  formulation 
within  the  subroutine  READPV.  The  triple  DO  loop  that  retrieves  the  element  numbers 
along  the  nonreflecting  boundary  involves  the  majority  of  the  additional  CPU  time. 
The  extra  storage  requirements  for  the  nonreflecting  boundary  are  minimal.  For  a 
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given  number  of  nodes,  N,  on  a  nonreflecting  boundary,  the  number  of  extra  variables 
stored  (the  size  of  the  array  BIMP)  in  the  program  would  be  2(r(  -  1). 

For  all  of  the  tests  conducted  in  this  study,  the  use  of  the  nonreflecting  bound¬ 
ary  hcis  required  a  negligible  increase  in  the  CPU  time.  Comparing  the  solution  time 
for  a  small  mesh  with  a  nonreflecting  boundary  and  the  solution  time  for  a  smail  mesh 
with  fixed/free  boundaries,  the  increase  in  CPU  time  was  less  than  3%.  Sometimes 
the  nonreflecting  boundary  solution  actually  required  less  time  than  the  solution  with 
the  fixed/free  boundary.  A  time  comparison  between  the  nonreflecting  boundary  solu¬ 
tion  and  the  large  mesh  solution  would  be  inapplicable.  The  solution  time  is  roughly 
proportional  to  the  number  of  nodes  in  a  mesh.  For  example,  if  the  number  of  nodes 
in  the  large  mesh  can  be  reduced  to  half  the  number  using  a  nonreflecting  boundary, 
the  solution  time  will  also  be  reduced  by  one-half.  The  savings  in  the  CPU  time  would 
depend  solely  upon  the  original  dimensions  of  the  large  mesh.  In  other  words,  the 
larger  the  original  mesh,  the  greater  the  savings  in  solution  time  using  a  small  mesh 
with  a  nonreflecting  boundary. 

The  stability  of  the  solutions  was  not  adversely  affected  in  all  of  the  wave  prop¬ 
agation  tests.  However,  the  maximum  allowable  time  step  that  could  be  used  in  the 
tests  was  not  determined.  Thus,  the  nonreflecting  boundary  can  be  deduced  not  to 
affect  the  stability  limit  of  a  problem  to  an  appreciable  extent,  if  at  all.  The  same  time 
step  for  a  large  mesh  can  be  used  for  a  smrdl  mesh  with  a  nonreflecting  boundary. 

This  study  has  shown  that  the  nonreflecting  boundary  almost  perfectly  absorbs 
stress  waves  that  impinge  normal  to  the  boundary.  The  effectiveness  of  the  nonreflect¬ 
ing  boundary  does,  in  fact,  decrease  with  an  increase  in  the  angle  of  incidence  of  the 
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stress  wave.  Therefore,  for  a  problem  involving  two-dimensional  wave  propagation,  the 
mesh  should  be  designed  such  that  the  nonreflecting  boundary  is  oriented  normal,  if 
possible,  to  the  wave  source. 

Axisymmetric  problems  can  be  difficult  to  model  for  a  nonreflecting  boundary 
that  is  parallel  to  the  axis  of  symmetry.  For  the  one-dimensional  wave  propagation  tests 
employing  an  axisymmetric  bar,  satisfactory  results  were  obtained  when  the  nonreflect¬ 
ing  boundary  parallel  to  the  axis  of  symmetry  was  placed  at  a  distance  corresponding 
to  an  R/A  value  greater  than  2.0.  In  a  two-dimensional  axisymmetrical  problem,  the 
wavelengths  of  the  radiating  stress  waves  cannot  be  calculated  accurately,  especially 
in  an  inelastic  medium.  Accordingly,  the  value  of  R/A  cannot  be  determined.  Thus, 
the  minimum  distance,  or  radius,  the  nonreflecting  boundary  can  be  from  the  axis  of 
symmetry  is  unknown  for  the  two-dimensional  case. 

Material  properties  that  were  analyzed  in  the  one-dimensional  wave  propaga¬ 
tion  tests  included  the  effect  of  Poisson’s  ratio  and  an  anisotropic  medium  on  the 
effectiveness  of  the  nonreflecting  boundary.  It  was  concluded  that  the  effectiveness  of 
the  nonreflecting  boundary  significantly  decreases  as  Poisson’s  ratio  nears  its  maximum 
value  of  0.5.  The  nonreflecting  boundary  performs  very  well  for  values  of  Poisson’s  ra¬ 
tio  between  0  and  0.4.  The  nonreflecting  boundary  is  also  very  effective  for  a  simple 
problem  involving  an  anisotropic  medium.  However,  the  ability  of  the  boundary  to 
absorb  stress  waves  in  a  two-dimensional  anisotropic  medium  is  still  undetermined  due 
to  the  complexity  of  the  problem. 

The  material  laws  used  in  this  study  were  material  law  6,  a  perfectly  elastic 
medium,  and  material  law  9,  the  AFWL  “engineering"  model.  The  remaining  material 
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laws  were  not  included  in  the  study.  The  biajcial  elastic-pl^ulic  material  laws,  applicable 
only  for  the  continuum  elements,  are  similar  to  the  perfectly  elastic  material  law  at  low 
values  of  stress.  The  inelastic  portion  of  the  elastic-plastic  material  law  will  be  similar  to 
a  multi-segmented  AFWL  “engineering”  law  at  higher  levels  of  stress.  The  endochronic, 
viscopl£istic,  and  the  Weidlinger  Associate’s  cap  materiaJ  law,  which  presently  exist 
in  the  code,  are  not  recommended  or  verified.  Further  testing  of  the  nonreflecting 
boundary  should  be  conducted  including  these  unverified  material  laws. 

The  two-dimensional  wave  propagation  tests  were  performed  using  the  AFWL 
“engineering”  model,  exclusively.  The  large  magnitude  of  stress  applied  to  the  system 
resulted  in  a  very  nonlinear  response.  The  dilatational  impedence  and  shear  impedence 
for  the  tests  corresponded  to  the  first  loading  segment,  or  the  “elastic”  section,  of 
the  hydrostatic  pressure  versus  volumetric  strain  curve  for  the  material.  The  use  of 
the  elastic  dilatational  impedence  and  elastic  shear  impedence  for  the  nonrefiecting 
boundary  gave  the  best  overall  results.  The  values  of  the  impedences  are  directly 
related  to  the  stress  wave  speed  in  the  material.  If  better  approximations  for  the 
stress  wave  speeds  in  a  problem  can  be  determined,  the  impedences  can  be  adjusted 
accordingly  by  changing  the  values  of  Young’s  modulus  and  Poisson’s  ratio.  Although 
not  examined  in  this  study,  this  adjustment  should  increase  the  effectiveness  of  the 
nonrefiecting  boundary. 

The  nonrefiecting  boundary  is  only  applicable  for  problems  involving  continuum 
elements  along  a  boundary.  The  interior  of  the  mesh,  however,  can  consist  of  any  type 
of  element.  Of  all  the  continuum  elements  in  the  SAMSON2  finite  element  library. 
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the  4NQ  works  best  with  the  iionreflecting  boundary.  This  was  evident  in  the  one¬ 
dimensional  wave  propagation  tests.  The  results  from  the  one-dimensional  tests  using 
the  5NT  and  6NT,  although  not  contained  in  the  report,  were  the  least  accurate  of  the 
continuum  elements. 

The  nonreflecting  boundary  cannot  be  used  in  a  static  or  quasi-static  problem. 
This  implies  that  the  use  of  the  dynamic  relaxation  feature  of  SAMSON2  is  prohib¬ 
ited.  When  utilizing  the  nonreflecting  boundary,  the  value  of  the  low-frequency  mass 
proportional  damping  parameter  should  be  zero. 

The  nonreflecting  boundary  should  not  be  placed  adjacent  to  a  load  line.  The 
nonreflecting  boundary  is  only  capable  of  absorbing  impinging  stress  waves  and  not 
the  direct  stress  resulting  from  the  direct  application  of  load. 

The  lone  application  of  Rayleigh  waves,  without  the  presence  of  dilatational 
waves  and  shear  waves,  has  not  been  included  in  this  study.  There  are  two  reasons 
for  this  exclusion.  First,  the  displacement  time  history  that  describes  a  Rayleigh  wave 
solution  is  cumbersome  to  discretize  for  a  SAMSON2  input  file.  Second,  and  most  im¬ 
portantly,  is  that  for  the  two-dimensional  wave  propagation  tests,  Rayleigh  waves  were 
also  generated  during  the  analyses  in  addition  to  the  two  body  waves.  The  displace¬ 
ment  and  stress  distribution  in  the  mesh  will  be  affected,  in  part,  by  the  propagation 
of  these  Rayleigh  waves.  If  the  nonreflecting  boundary  solution  compares  well  to  the 
large  mesh  solution,  the  nonreflecting  boundary  can  be  assumed  to  be  capable  of  also 
absorbing  Rayleigh  waves. 

Fulfilling  the  second  objective  of  the  research,  the  pos.sible  errors  that  exist  in 


the  SAMSON2  code  have  been  found.  The  anomalies  discovered  in  the  SAMSON2 
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code  and  SAMS0N2  User’s  Manual  by  the  authors  are  located  in  Appendix  B  along 

with  recommended  corrections. 

The  following  list  provides  a  summary  of  the  recommendations  on  the  use  of 

the  nonreflecting  boundary  presented  in  this  chapter. 

1.  The  nonrefiecting  boundary  should  be  oriented  normal  to  the  wave  source  if 
possible. 

2.  The  nonreflecting  boundary  parallel  to  the  axis  of  symmetry  should  be  located 
such  that  the  distance  from  the  boundary  to  the  cixis  corresponds  to  an  R/A  ratio 
greater  than  2.0  if  the  wavelength  of  the  stress  wave  can  be  approximated. 

3.  When  the  AFWL  “engineering”  material  model  is  employed  with  a  multi-segmented 
loading  curve,  the  Young’s  modulus  and  Poisson’s  ratio  should  correspond  to  the 
first  loading  segment  if  a  typical  wave  speed  for  the  problem  cannot  be  accurately 
determined. 

4.  Use  of  the  4NQ  continuum  element  yields  the  best  results  for  the  nonreflecting 
boundary. 

5.  The  nonreflecting  boundary  is  not  applicable  for  a  static  problem  (i.e.,  the  use 
of  dynamic  relaxation  is  not  allowed). 

6.  The  nonreflecting  boundary  should  not  be  placed  adjacent  to  a  load  line. 
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CARD  lOA 

LOAD  LINE  CONTROL  READPV 

FORMAT(7I5,5X,E10.3) 

Columns 

Variable 

Description 

1-5 

I 

Load  line  number 

6-10 

NDNOD 

Number  of  nodes  on  load  line 

11-15 

IVOL(l) 

=  -2,  pressure  line  for  axisymmetric  problems 
=  -1,  pressure  line  for  plane  problems 
=  0,  initied  velocity  line 
=  1,  force  line 

=  2,  displacement  history  line 
=  3,  axisymmetric  nonreflecting  boundary 
=  4,  plane  nonreflecting  boundary 

16-20 

rVOL(2) 

Used  only  when  IVOL(l)  =  0 
=  0,  initial  velocities  are  normal  and  tangent  to 
the  load  line 

=  1,  initial  velocities  are  in  the  global  x-  and 
y-componcnts 

21-25 

IVOL(3) 

Force  or  displacement  component  code.  Used 
only  when  IVOL(l)  =  1  or  2 
=  1,  x-component  is  specified 
=  2,  y-component  is  specified 

26-30 

INTI 

Node  generation  increment 

31-35 

NPT 

Number  of  points  used  to  define  the  piecewise 
linear  pressure,  force,  or  displacement  history. 

Not  used  if  IVOL(l)  =  0,  3,  or  4 

36-40 

Blank 

41-50 

Thick 

Thickness  for  planar  problem  pressure 

#  modifications  are  printed  in  bold  type 
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Appendix  B 

ANOMALIES  IN  THE  SAMSON2  CODE  AND  USER’S  MANUAL 

The  latest  version  of  the  SAMSON2  code  was  received  by  the  authors  late  in 
the  period  of  the  research.  All  of  the  errors  that  were  discovered  in  the  outdated 
SAMSON2  code  -  which  was  used  primarily  in  the  research  -  were  inapplicable  upon 
receiving  the  new  version  of  the  code.  The  only  anomaiy  that  was  found  in  the  latest 
version  of  the  code  is  located  in  subroutine  READMA.  The  last  executable  line  in  the 
subroutine  before  the  RETURN  statement  computes  the  value  of  LADD.  The  last  term 
in  the  line  subtracts  one  from  the  sum  of  INDXE(NUMMAT)  plus  LE2  plus  LADD2. 
To  be  precise,  two  should  be  subtracted  from  the  sum  instead  of  one.  However,  the 
only  advantage  to  this  correction  is  a  negligible  savings  in  the  storage  requirement  for 
the  program.  The  correction  is  simply  cosmetic. 

There  were  typographic  mistakes  in  the  input  instructions  section  of  the  revised 
SAMSON2  User’s  Manual  [13].  Several  card  format  statements  in  the  User’s  Manual 
either  do  not  correlate  with  the  card  column  specifications  or  do  not  match  the  format 
statements  in  the  SAMSON2  code.  Some  of  these  anomrdies  in  the  format  statements 
are  purely  academic;  they  do  not  affect  the  column  specifications  at  all  or  are  simply 
a  variation  in  the  type  of  variable  (Ex.  FIO.O  in  the  code  might  be  listed  as  FlO.3  in 
the  User’s  Manual).  Those  anomalies  that  do  affect  the  column  specifications  should 
be  corrected  in  the  next  revision  of  the  SAMSON2  User’s  Manual  and  are  discussed  in 


this  appendix. 
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The  first  anomaly  is  on  card  4  which  contains  the  program  control  data.  The 
last  variable,  IPPLF,  is  not  read  in  the  execution  of  the  code.  The  reason  is  unknown  to 
the  authors.  The  column  specification  for  the  alphanumeric  title  in  card  9B,  “motion 
control  output”,  should  be  16-31  instead  of  11-23.  Likewise  for  card  9C,  “element 
variable  output”,  the  column  specification  for  the  title  should  be  16-31  instead  of  16- 
32.  The  most  significant  error  found  in  the  User’s  Manual  was  found  on  card  9D, 
“picture  output  request”.  The  format  statement  should  read  215  instead  of  2110.  Also, 
the  column  specifications  for  the  variables  NTSTEP  and  K  should  be  1-5  and  6-10, 
respectively. 

The  input  cards  that  do  not  have  errors  in  their  format  statements  are  the 
following:  card  1,  “title”;  card  3,  “master  control  data”;  card  4,  “program  control 
data”;  card  5C,  “material  law  parameter  data”;  card  8,  “nodal  boundary  condition 
data”;  card  9B,  “motion  component  output”;  card  9C,  “element  variable  output”; 
card  9F,  “node  point  output”;  card  9G,  “element  output”;  card  lOB-1,  “load  line 
history”;  card  IOC,  “load  line  nodes”;  and  all  of  the  slideline  cards.  The  correct  format 
statements  for  the  cards  not  mentioned  above  can  be  found  in  the  latest  version  of  the 


SAMSON2  code  used  in  this  study. 


Appendix  C 

MODIFIED  SUBROUTINES  IN  THE  SAMSON2  CODE 
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The  subroutines  that  have  been  modified  by  the  authors  are  listed  in  this 
appendix.  The  following  FORTRAN  listings  have  been  changed  from  the  AFWL’s 
HP9000  to  be  compatible  to  the  IBM  3090  computer.  The  only  differences  occur  in 
the  initial  lines  at  the  beginning  of  each  subroutine.  The  common  block,  character, 
equivalence  and  dimension  statements  in  the  IBM  3090  version  are  substituted  for  the 
include  statements  in  the  HP9000  version  of  the  code.  The  modified  executable  lines 
are  printed  in  bold  type. 


r>  n  no 


SUBROUTINE  FREEF2  (XC.YC.Xl.V.FORCD.STRS.lNDXS.IP.NPT.NDNOD. 

1  IVOL.KPRES.PT.THICK,  BIMP) 

C - - - 

C 

C  CALCULATES  VELOCITY  AND  PRESSURE  BOUNDARY  CONDITIONS 
C 

C - 

COMMON/DYNAM/  DELT.TIME.MXSTEP.NTSTEP.NTSTPl.OMEGA.WMAX.CWAVE, 
-  CBULK,CSHEAR,C1DAMP,TIMEMX,DTMIN.DTMAX,MXSUB,P0ISR 

COMMON/RTIME/  DATE.  RUNTIM,  CREATR,  VRSION 
CHARACTER*8  DATE,  RUNTIM,  CREATR,  VRSION 
C 

COMMON/NUM/  NUMNP.NUMEL,NUMMAT.NUMDIS.NOGREE,NPRES,NSLID,MEQ, 

*  DELTMX.LE2,NPLOT 

COMMON/KNTRL/ISYMF,MITSF,IUN5FITllF,IMSHF,IHISF, 

+  IQARF.IRSTF 

DIMENSION  KONTRL(IO) 

EQUIVALENCE  (KONTRL(l).ISYMF) 

C 

DIMENSION  XC(1).  YC(1),  Xl(l),  V(l).  FORCD(l).STRS(l). 

1  INDXS(l),  IVOL(5),  KPRES(l),  PT(1),  BIMP(l) 

C 

IF(IVOL(l)  .EQ.  0)  RETURN 
NDN  =  NDNOD-l 
TOTFX  =  0.0 
TOTFY  =  0.0 
DVOL  =  0.0 
C 

C****  COMPUTE  NODAL  FORCES  FOR  PRESSURE  LOADINGS  ***• 

C»***  GET  OLD  VOLUME  FOR  PRESSURE  LINE  OR  OLD  TOTAL  DISPLACEMENT  ***• 
C****  GET  OLD  PRESSURE  OR  FORCE  **** 

C****  GET  OLD  EXTERNAL  WORK  **** 

C****  GET  OLD  IMPULSE  *•** 

C 

K  =  INDXS(NUMEL+2)+IP-l 
VOLD  =  STRS(K) 

C 

L  =  INDXS(NUMEL+3)+IP-l 
POLD  =  STRS(L) 

C 

M  =  INDXS(NUMEL+5)+IP-l 
WEXT  =  STRS(M) 

C 

MM  =  INDXS(NUMEL+4)+IP-l 
WIM  =  STRS(MM) 

***•  PRESSURE.FORCE,  DISPLACEMENT  HISTORY 
CALL  PRESS(  NPT,  PT,  P,  TIME) 

IF  (IVOL(l)  .EQ.  1)THEN 

***•  NODAL  FORCE  LOADING  ***• 
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C - 

J  =  IVOL(3) 

C 

DO  510  II  =  l.NDNOD 
INI  =  KPRES(ll) 

IN2  =  INl*2-2-t-J 
FORCD(IN2)  =  FORCD(IN2)+P 
C 

C****  COMPUTE  TOTAL  FORCE  *  DELT  X 
TOTFX  =  TOTFX  +  V(IN2) 

510  CONTINUE 
C 

C****  STORE  NEW  FORCE, IMPULSE,  TOTAL  WORK  IN  STRS  ARRAY 
C 

STRS(L)  =  P 

STRS(MM)  =  WIM+(P+POLO)*NDNOD*.5*DELT 
STRS(M)  =  WEXT+(P  +  POLD)*.5*TOTFX*OELT 
C 

ELSEIF  (IVOL(l)  .EQ.  2)  THEN 

C - 

C***’  DISPLACEMENT  TIME  GIVEN 
C‘***  COMPUTE  VELOCITY  •*** 

C - 

J  =  IVOL(3) 

DO  610  II  =  l,NDNOD 
INI  =  KPRES(ll) 

IN2  =  INl*2-2+J 
V(IN2)  =  (P-X1(IN2))/DELT 
Xl(IN2)  =  P 
610  CONTINUE 
C 

ELSEIF  ((IVOL(l)  .EQ.  -1)  .OR.  (IVOL(l)  .EQ.  -2))  THEN 
C - 

C»***  pressure  line 

c - 

DO  100  II  =  1,NDN 
C 

C****  FIND  NEXT  TWO  NODES  NO  LOAD  LINE 
INI  =  KPRES(ll) 

IN2  =  KPRES(ll  +  l) 

XCl  =  XC(INl)+Xl(2*INl-l) 

XC2  =  XC(IN2)+X1(2*IN2-1) 

YCl  =  YC(IN1)+X1{2*INI) 

YC2  =  YC(IN2)+X1(2*IN2) 

SS  =  XCI-XC2 
CC  =  YC2-YC1 
C 

C****  COMPUTE  CHANGE  IN  VOLUME  **** 

PP  =  (YC(IN2)-YC(IN1))*(X1(2*INI-1)+X1(2*IN2-1)) 

1  -(XC(IN2>XC(IN1))*(X1(2*INI)+XI(2*IN2)) 

2  +X1(2*IN1-1)*X1(2*IN2)-XI(2*IN1)*X1{2*IN2-1) 

FX  =  0.5*CC*P*THICK 


n  n  n  n  rt 
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FY  =  0.5*SS*P*THICK 
IF  (IVOL(l)  .IME.  -1)  THEN 
C 

C*‘**  AXISYMMETRIC  LOADING 
CC  =  0.5*(XC1+XC2) 

PP  =  PP*CC 
FX  =  FX*CC 
FY  =  FY*CC 
ENDIF 

OVOL  =  DVOL-0.5*PP 
C 

C**»‘  STORE  EXTERNAL  FORCE  ***• 

TOTFX  =  TOTFX  +  FX 
TOTFY  =  TOTFY+FY 
FORCD(2*INl-l)  =  FORCD(2*INl-l)+FX 
FORCD(2*IN2-l)  =  FORCD(2*IN2-l)+FX 
FORCD(2»INl)  =  F0RCD(2*IN1)  +  FY 
F0RCD(2*IN2)  =  FORCD(2*IN2)+FY 
100  CONTINUE 
C 

C****  STORING  VOLUME  CHANGE, WORKDONE  BY  PRESSURE,  PRESSURE,  IMPULSE 
C»***  IN  STRS 
C 

STRS(K)  =  DVOL 

STRS(M)  =  WEXT+0.5*(POLD  +  P)*(DVOL-VOLD) 

STRS(L)  =  P 
IF(P  .NE.  0.0)  THEN 

TOTFX  =  SQRT(TOTFX*TOTFX+TOTFY*TOTFY)*DELT*(1.0+POLD/P)V5 

/• 

V. 

C****  LAST  TERM  IS  AN  APPROXIMATION  TO  F+FOLD  **** 

STRS(MM)  =  WIM  +  TOTFX 
ENDIF 
C 

ELSE 


NONREFLECTING  BOUNDARY 


DO  200  II  =  l.NON 

FIND  NEXT  TWO  NODES  NO  LOAD  LINE 
INI  =  KPRES(ll) 

IN2  =  KPRES(I1+1) 

XCl  =  XC(IN1)+X1(2*INM) 

XC2  =  XC(IN2)+X1(2*IN2-1) 

YCl  =  YC(IN1)-FX1(2*IN1) 

YC2  =  YC(IN2)+X1(2*IN2) 

SS  =  XC2-XC1 
CC  =  YC2  YCl 

c 

REFP=BIMP(NDNOD+2*ll-l) 

REFV=BIMP(NDNOD+2*ll) 

C 
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COMPUTE  CHANGE  IN  VOLUME 

PP  =  (YC(IN2)-YC(IN1))*(X1(2*INI-1)+X1(2*IN2-1)) 

1  -(XC(IN2)-XC{IN1))*(X1(2*INI)+X1(2*IN2)) 

2  +Xl(2*INl-l)*Xl(2*IN2^Xl(2*tNl)*Xl(2*IN2-l) 

C 

V1X=V(2*INM) 

V1Y=V{2*IN1) 

V2X=V(2*IN2-1) 

V2Y=V(2*IN2) 

C 

FX1=-(2./6*VIX  +1./6.*V2X)*(CC*REFP  +SS*REFV)*THICK 
FX2=-(1./6.*V1X  +2./6.*V2X)*(CC*REFP  +SS*REFV)*THICK 
FYI=-(2./6.*VlY  +1./6.*V2Y)*(SS*REFP  +CC*REFV)*THICK 
FY2=-(1./6.*VIY  +2./6.*V2Y)*(SS*REFP  +CC*REFV)*THICK 
IF  (IVOL(l)  .EQ.  3)  THEN 
C 

C****  AXISYMMETRIC  LOADING  **** 

CC  =  0.5*(XC1+XC2) 

PP  =  PP*CC 
FXl  =  FX1*CC 
FYl  =  FY1*CC 
FX2  =  FX2*CC 
FY2  =  FY2*CC 
ENDIF 

DVOL  =  DVOL+0.5*PP 
C 

C****  STORE  EXTERNAL  FORCE  *♦** 

TOTFX  =  TOTFX+FXH-FX2 
TOTFY  =  TOTFY+FY1+FY2 
fOilCO(2*INll)  =  FORCO(2*INl-l)+FXl 
FORCO(2*IN2-l)  =  F0RC0(2*IN2.I)+FX2 
FORCO(2*IN1)  =  FORCO(2*INl)+FYl 
FORCO(2*IN2)  =  FORCD(2*IN2)+FY2 
C 

WEXT  =  WEXT  +  ABS(FX1*XI(2*IN1-1)) 

WEXT  =  WEXT  +  ABS(FX2*XI(2*IN1)) 

WEXT  =  WEXT  +  ABS(FY1*X1(2*IN2-1)) 

WEXT  =  WEXT  +  ABS(FY2*Xl(2*IN2)) 

C 

200  CONTINUE 
C 

C«*««  STORING  VOLUME  CHANGE.WORKDONE  BY  PRESSURE.  PRESSURE,  IMPULSE 
C****  IN  STRS 
C 

STRS(K)  =  DVOL 
STRS(M)  =  WEXT 

TOTFX  =  SQRT(TOTFX*TOTFX+TOTFY*TOTFY)*DELT 
C 

C****  LAST  TERM  IS  AN  APPROXIMATION  TO  F+FOLO 
STRS(MM)  =  WIM  +  TOTFX 

C 

ENDIF 


RETURN 

END 


o  o 
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SUBROUTINE  FREEFD  (  XC.  YC.  XI,  V.  FORCD, 

1  SIRS,  INDXS.  IP.  NPT,  NDNOD, 

2  IVOL,  KPRES,  PT,  THICK,  BIMP) 


C***»  THIS  ROUTINE  WILL  COMPUTE  THE  LOADING  CONDITIONS  AS  A  FUNCTION**** 
C****  OF  TIME.  INITIAL  CONDITIONS  ARE  COMPUTED  IF  FREEFD  IS  CALLED  **** 
C*.**  pressure  LOADINGS  ARE  COMPUTED  IF  FREEF2  IS  CALLED. 

^****  ***« 

C  1.  COMPUTES  NODAL  FORCES  FOR  A  GIVEN  PRESSURE  LINE  WITH 
C  PIECEWISE  LINEAR  PRESSURE  TIME  HISTORY 

C  2.  COMPUTES  NODAL  FORCES  THAT  ARE  DIRECTLY  SPECIFIED  WITH 
C  PIECEWISE  LINEAR  TIME  HISTORY 

C  3.  SETS  NONZERO  INITIAL  VELOCITIES  FOR  IMPULSE  LOADING 
C  NORMAL  TO  LOAD  LINE  OR  IN  GLOBAL  X,Y  DIRECTIONS 

C  4.  COMPUTES  NONZERO  PRESCRIBED  DISPLACEMENTS  FROM  PIECEWISE 
C  LINEAR  INPUT 

C 

C  XC.YC  =  GLOBAL  (ORIGINAL)  X.Y  COORDINATE  ARRAYS 
C  BIMP  =  BOUNDARY  IMPEOENCE  ARRAY 
C  XI  =  nodal  DISPLACEMENT  ARRAY 
C  V  =  NODAL  VELOCITY  ARRAY 
C  FORCD  =  GLOBAL  EXTERNAL  FORCE  ARRAY  .  F-EXT 
C  STRS  =  ELEM.STORAGE  ARRAY 

C  INDXS  =  LOCATOR  FOR  STRS 

C  IP  =  LOAD  LINE  NO. 

C  NPT  =  NO. OF  PTS.FOR  PIECEWISE  LINEAR  DESCRIPTION  OF  FORCES 
C  OR  DISPLACEMENTS  ON  LOAD  LINE 

C  NONOO  =  NO.OF  NODES  ON  LOAD  LINE 
C  PT  =  VALUES  OF  FUNCTION  AND  TIME  AT  NPT  POINTS 
C  IVOL  =  CONTROL  FOR  LOAD  LINE  *•  REFER  TO  READPV 
C  KPRES  =  NODES  DEFINING  THE  LOAD  LINE 
C 

C****  NUMDIS  NUMBER  OF  NODES  WITH  DISPLACEMENT  PRESSURE 
C***.  PRESCRIBED 

C***«  NNODE  NUMBER  OF  NODES  IN  PROBLEM 
C****  XC(N)  X  COORDINATE  OF  NODE  NUMBER  N 
C****  YC(N)  Y  COORDINATE  OF  node  number  N 
C*«*«  NDNOD  NUMBER  OF  NODES  ON  LOAD  LINE 
C****  IVOL(l)  =  -IPRESSURE  TIME  INPUT  AXISYMMETRIC 
C****  IVOL(l)  =  0  PRESSURE  TIME  INPUT  PLANE 
C****  IVOL(l)  =  1  INITIAL  VELOCITY  INPUT 

C****  IVOL(2)  =  0  LOAD  APPLIED  NORMAL  AND  TANGENT  TO  LOAD  LINE 

C****  IVOL(2)  =  1  LOAD  APPLIED  IN  GLOBAL  X  AND  Y  DIRECTIONS 

C***.  IV0L(3)  USED  FOR  FORCE  LOADING 

C****  =  0  FORCE  IS  APPLIED  IN  GLOBAL  X  DIRECTION 

C****  =  I  FORCE  IS  APPLIED  IN  GLOBAL  Y  DIRECTION 

C****  STRS  ELEMENT  STRESS  ARRAY 

C****  INDXS  INDEX  TO  ELEMENT  STRESS  ARRAY  (STRS) 

C - 

C 
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COMMON/DYNAM/  DELT.TIME.MXSTEP.NTSTEP.NTSTPl, OMEGA, WMAX.CWAVE, 
+  CBULK,CSHEAR,C1DAMP,TIMEMX,DTMIN,DTMAX.MXSUB,P0ISR 

COMMON/RTIME/  DATE,  RUNTIM,  CREATR,  VRSION 
CHARACTER'S  DATE,  RUNTIM,  CREATR,  VRSION 
C 

COMMON/NUM/  NUMNP,NUMEL,NUMMAT,NUMDIS,NDGREE,NPRES,NSL1D,MEQ, 

•  DELTMX,LE2,NPLOT 

COMMON/KNTRL/ISYMF,MITSF,IUN5F,lTllF,IMSHF,IHISF, 

+  IQARF,IRSTF 

DIMENSION  KONTRL(IO) 

EQUIVALENCE  (KONTRL(l),ISYMF) 

C 

DIMENSION  XC(1),  YC(1),  Xl(l),  V(l),  FORCD(l),  STRS(l), 

1  INDXS(l),  IVOL(S),  KPRES(l),  PT(1),  BIMP(l) 

C 

C****  CONTOL  FOR  INITIAL  VALUE  INPUT  •*** 

IF  (IVOL(l)  .EQ,  0)  THEN 
C 

C****  INITIAL  VELOCITY  INPUT  ***♦ 

IF(IVOL(2)  .EQ.  1)  THEN 
C 

C****  GLOBAL  VELOCITY  LOADING  ♦*** 

DO  176  II  =  l,NDNOD 
INI  =  KPRES(ll) 

IN2  =  IN1*2 

V(IN2-1)  =  V(IN2-1)+PT(1) 

V(IN2)  =  V(IN2)+PT(2) 

176  CONTINUE 
C 

ELSE 

NON  =  NDNOD-1 
DO  190  II  =  1,NDN 
C 

C*’**  FIND  CONSECUTIVE  NODES  ON  LOAD  LINE  **** 

INI  =  KPRES(ll) 

IN2  =  KPRES(11  +  1) 

CC  =  XC(IN2)-XC(IN1) 

SS  =  YC(IN2)-YC(IN1) 

AL  =  SQRT(SS*SS+CC*CC) 

CC  =  CC/AL 
SS  =  SS/AL 
C 

C****  APPLY  HALF  VELOCITY  AT  EACH  NODE  •*** 

VT  =  0.5*PT(1) 

VN  =  0.5*PT(2)  ^ 

VX  =  VT*CC+VN*SS 
VY  =  VT*SS-VN*CC 
C 

IF  (ll.EQ.l  .AND.  PT(3).NE.0.0)  THEN 
C 

C****  IF  FIRST  NODE  APPLY  FIRST  NODE  FACTOR  **♦* 

V(2*IN1-1)  =  V(2*IN1-1)  -I-  VX*PT(3) 


126 


V(2*1N1)  =  V(2*1NI)  +  VY*PT(3) 

V(2*IN2-1)  =  V(2*1N2-1)  +  VX 
V(2*IN2)  =  V(2*IN2)  +  VY 
ELSEIF(I1.EQ.NDN  .AND.  PT(4).NE.0.0)  THEN 
C 

c»***  apply  last  node  factor  *’•* 

V(2*IN1-1)  =  V(2*1N1-1)  +  VX 
V(2*IN1)  =  V(2*IN1)  +  VY 

V(2*IN2-1)  =  V(2*IN2-1)  +  VX*PT(4) 

V(2*IN2)  =  V(2*IN2)  +  VY*PT(4) 

ELSE 

V^2•|Nl-l)  =  V(2*IN1-1)  +  VX 
V(2*IN1)  =  V(2*IN1)  +  VY 
V(2*IN2-1)  =  V{2*IN2-1)  +  VX 
V(2*IN2)  =  V(2*IN2)  +  VY 
ENDIF 
C 

190  CONTINUE 
ENDIF 
ENDIF 
C 

C**..  pressure  TIME  LOADING 

CALL  FREEF2  (XC.YC,Xl.V,FORCD,STRS.INDXS.IP.NPT,NDNOD,IVOL, 
1  KPRES.PT, THICK.  BIMP  ) 


C 


RETURN 

END 
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SUBROUTINE  PREPRC 

C - 

C  1.  CONTROLS  READING  OF  INPUT  DATA  FILE 

C  2.  SETS  POINTERS  (LVARIA8LES)  FOR  ALL  ARRAYS  AND  CHECKS  THAT 
C  Q-ARRAY  SIZE  IS  CONSISTENT  WITH  AVAILABLE  CORE. 

c 

c. - - - - - 

DIMENSION  IQ(1) 

EQUIVALENCE  (Q{1),IQ(1)) 

COMMON  /MAXQ/  MAXQ 
COMMON  /ARRAY/  Q(l) 

C 

COMMON/DYNAM/  DELT.TIME.MXSTEP.NTSTEP.NTSTPl.OMEGA.WMAX.CWAVE. 

+  cbulk,cshear,cidamp.timemx,dtmin.dtmax.mxsub,poisr 

COMMON/RTIME/  DATE,  RUNTIM,  CREATR,  VRSION 
CHARACTER’S  DATE,  RUNTIM,  CREATR,  VRSION 
C 

C  WARNING -  DO  NOT  CHANGE  LOCAT  EXCEPT  AFTER  LTOTAL  OR 

C  SUBROUTINE  DEBUG  WILL  NOT  WORK 

COMMON/LOCAT/  LFINT  ,LXC  ,LYC  ,LNODDI,LIX  ,LXl  ,LVl  , 

+  LA1  ,lFORCD,LXO  ,LDTN0D,LSMASS,LMAP  .LICTYP.LDTGRP.LANGLE, 
+LINDXE,LE  ,LIPLT  ,LELOUT,LOUT  .LNP0UT,LINDXP,LPRS  .LINDXL, 

+  LSLO  ,LINDXS,LSTRS  ,LPSU  ,LT  ,LUU  ,LIXX  ,LIYY  , LTOTAL 
DIMENSION  LPTR(35) 

EQUIVALENCE  (LPTR, LFINT) 

C 

COMMON/OUTPA/  NPRU.NPRS,NPFREQ,NPIC,NKN,MXN,NPRUS 
COMMON/OUTPC/  TITLE 
CHARACTER  TITLE’80 
C 

COMMON/NUM/  NUMNP,NUMEL,NUMMAT,NUMDIS,NDGREE,NPRES.NSLID,MEQ. 

♦  DELTMX,LE2.NPLOT 

COMMON/KNTRL/ISYMF.MITSF,IUN5F,ITllF,IMSHF,IHISF, 

+  IQARF,IRSTF 

DIMENSION  KONTRL(IO) 

EQUIVALENCE  (KONTRL(l),ISYMF) 

C 

COMMON/UNITS/LTNM,LLNM,LFNM,LWNM,GMAG.LPNM,LVNM.LDNM 
COMMON/UNITC/  TIMENM,LENGNM,FORCNM,WGTNM,PRESNM  ,VELNM, 

+  DENSNM 

CHARACTER*4  TIMENM,LENGNM,FORCNM,WGTNM,PRESNM*12,VELNM*9, 

+  DENSNM*20 
CHARACTER  COMMNT’80 
C 

C”**  PRINT  BANNER,DATE,  TIME  AND  VRSION 
C 

CREATR  =  ’SAMSON2’ 

VRSION  =  ’0.8’ 

CALL  RTIME(DATE, RUNTIM) 


onn  rvr>o  on  nr>r> 
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WRITE  (6,1000) 

WRITE  (6,1001) 

WRITE  (6,1002)  CREATR,VRSION,DATE,  RUNTIM 
C 

C****  INITIALIZE 
C 

REWIND  1 
REWIND  2 
NDGREE  =  2 
NTSTEP  =  0 
NTSTPl = 1 
TIME  =  0.0 
DELTMX  =  9.99E10 
NKN  =  0 
C 

C****  READ/WRITE  USER  COMMENT  BLOCK.  TERMINATE  WITH  BLANK  STRING 
C 

DO  90  I  =  1,  100 
READ  (5.110)  COMMNT 
(F(COMMNT(l:2)  .NE.  '•/■)GO  TO  lOO 
IF(I  .EQ.  1)  WRITE(6.109) 

WRITE(6,111)  COMMNT(3;80) 

90  CONTINUE 
C 

100  BACKSPACE  5 
READ(5,110)  TITLE 

•'**  CHECK  FOR  RESTART  RUN  IN  FIRST  PART  OF  PROBLEM  TITLE  *  *  * 

if(T1TLE(l;7).EQ. ’RESTART’  .OR.  TITLE(l:7).EQ. ’RESTART’  )  THEN 
CALL  RESTAR(2) 

IF(IHISF.NE.O)  CALL  PSTAPE 

IF  THE  USER  HAS  REQUESTED  OUTPUT  TO  THE  COMMON  DATA  BASE  (ITllF 
IS  NONZERO)  THEN  WRITE  ALL  DISCRIPTIVE  INFORMATION  TO  TAPEll 

IF  (ITllF  .NE.  0)  THEN 

CALL  INIT11(Q(LIX),Q(LXC),Q(LMAP).Q(LIPLT)) 

END  IF 

PRINT  OUT  SELECTED  NODES  AND  ELEMENT  DATA  ON  RESTART  CYCLE 

CALL  OUTSNP(Q(LXC),Q(LYC),Q(LXl),Q(LVl),Q(LAl),Q(LINDXS), 

+  Q(LINDXE),Q(LSTRS),Q(LIX).Q(LE).Q(LMAP),Q(LICTYP)) 

NEXT  COMPUTATIONAL  CYCLE  IS  RESTART  CYCLE  +  1 

NTSTPl  =  NTSTEP  + I 
RETURN 
ENDIF 


WRITE  (6,120)  TITLE 


non 


C 

C****  READ  IIM  G-MAGNITUDE  AND  TIME.LENGTH, FORCE  UNITS  **♦* 

C 

CALL  READUN 
C 

C****  READ(WRITE)  PROBLEM  SIZE  PARAMETERS  ***** 

C 

READ  (5.130)  NUMNP.NUMEL.NUMMAT.NUMDIS.NPRES.NSLID.MXSTEP.DELT, 

1  TIMEMX.CIDAMP 

WRITE  (6,5)  NDGREE.NUMNP.NUMEL.NUMMAT.NUMDIS.NPRES.NSLID.MXSTEP, 
DELT.TIMENM.TIMEMX.TIMENM.ClDAMP 

•*•*  READ  CONTROL  PARAMETERS  ***** 

CALL  READKO 
C 

C****  SET  NSIZE  ARRAY  ***** 

C 

CALL  SETNSZ 
C 

C****  SET  ADDRESSES  ****** 

C 

MEQ  =  NDGREE*NUMNP 
LFINT  =  I 

LXC  =  LFINT  +  MEQ 
LYC  =  LXC  +  NUMNP 
LNODDI  =  LYC  +  NUMNP 
LIX  =  LNODDI  +  3*(NUMDIS+1) 

LXl  =  LIX  +  10*NUMEL 

LVl  =1X1  4- MEQ 

LAI  =  LVl  +  MEQ 

LFORCD  =  LAI  +  MEQ 
LXO  =  LFORCD  +  MEQ 
LOTNOD  =  LXO  +  MEQ 
IF(MITSF  .LE.  0)  THEN 
LSMASS  =  LDTNOD 
ELSE 

LSMASS  =  LDTNOD  +  NUMNP 
ENDIF 

LMAP  =  LSMASS  +  MEQ 
LICTYP  =  LMAP  +  NUMMAT 
LDTGRP  =  LICTYP  +  NUMMAT 
LANGLE  =  LDTGRP  +  NUMMAT 
LINDXE  =  LANGLE  +  NUMDIS  +  1 
LE  =  LINDXE  +  NUMMAT 
C 

C  MEMADJ  IS  A  CFTLIB  ROUTINE  TO  ADD  MEMORY  TO  THE  END  OF  BLANK 
C  COMMON  AT  RUN  TIME,  IT  RETURNS  lERR  =  0  IF  SUCCESSFUL 
C  ESTIMATE  AMOUNT  OF  CORE  NEEDED  AND  ADJUST  ACCORDINGLY 
C 

MORMEM  =  LE  +  NUMMAT*100  +  NPRES*400  +  NSLID*400  +  7*NUMEL 
+  +  7*(NPRS+NPRU)  +  2*NPIC 


r>  o  r>  n  r>  r> 


C  CALL  DYNAMEM(MORMEM,MAXQ) 

C 

CALL  ZERO  (Q.l.MAXQ) 

C 

C****  READ  IN  MATERIAL  DATA(LADD  IS  LENGTH  OF  E-ARRAY(El  AND  E2)) 
€**••  NODAL  COORDINATES  **** 

C****  ELEMENT  DATA 

C*’**  INITIALIZE  MULT.  INTEG.  TIMESTEP  DATA 

C****  BOUNDRY  CONDITIONS  •*** 

C****  OUTPUT  CODE  CARDS  **•* 

C 

CALL  READMA  (Q(LE),Q(LINDXE).Q(LMAP).Q(LICTYP).Q(LDTGRP),LADD) 
CALL  READNO  (Q(LXC),Q(LYC)) 

CALL  READEL  (Q(LIX).Q(LMAP),Q(LDTNOD).Q(LDTGRP)) 

CALL  INMITS  (Q(LIX).Q(LDTNOD).Q{LDTGRP)) 

CALL  READBC  (Q(LNODDI),Q(LANGLE)) 

C 

READ  (5,6)  NPFREQ.NPRU.NPRS.NPIC 
C 

LIPLT  =  LE  +  LADD 
LELOUT  =  LIPLT  +  NUMEL 
LOUT  =  LELOUT  +  6*NUMEL 
LNPOUT  =  LOUT  +  7*(NPRU  +  NPRS) 

C 

CALL  READOU  (Q(LOUT),Q(LNPOUT)) 

CALL  READSN 

IF(IHISF  .NE.  0)CALL  INITSN(Q(LXC),Q(LYC).Q{LE).Q(LMAP), 
Q(LINDXE),Q(LIX)) 

•***  AEAD  IN  LOADING  DATA  •*** 

LINDXP  =  LNPOUT  +  2*NPIC 
LPRS  =  LINDXP  +  NPRES 
LADD  =  0 

IF  (NPRES.NE.O)  THEN 
DO  1  I  =  1, NPRES 
L  =  LPRS  +  LADD 
IQ(LINDXP+I-1)  =  1  +  LADD 

CALL  READPV(Q(L),Q(L+5),Q(L+6),Q{L+6),LI.  Q(LE).Q(LINDXE), 
Q(LIX).Q(l+6)  ) 

LADD  =  LADD  +  LI 
1  CONTINUE 
ENDIF 

READ  IN  SLIDE  LINE  DATA 

LINDXL  =  LPRS  +  LADD 
LSLD  =  LINDXL  +  NSLID 
LADD  =  0 
IF  (NSLID.NE.O)  THEN 
DO  3  1  =  1, NSLID 
L  =  LSLD  +  LADD 


n  n  r>  nnno  rtrtnnrtrtnn  no 
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IQ(LINDXL+I-1)  =  1  +  LADD 
CALL  READSL(Q(L),Q(L+11).L1) 

LADD  =  LADD  +  LI 
3  CONTINUE 
ENDIF 

LINDXS  =  LSLD  +  LADD 
LSTRS  =  LINDXS  +  NUMEL  -t-  6 
C 

C****  MORE  BLANK  COMMON  WILL  BE  REQUIRED  AT  THIS  POINT,  CALL  DYNAMEM 
estimate  amount  and  ADJUST  ACCORDINGLY 
C**»*  LADD  -  SIZE  PRINTER  PLOT  HISTORY  ARRAYS 

C****  LADDS  -  AMOUNT  OF  STORAGE  IN  Q-ARRAY  EXCLUSIVE  OF  STRS  ARRAY 
C****  LQUES  -  QUESS  AT  FINAL  SIZE  OF  Q-ARRAY  INCLUDING  STRS  ARRAY 
*•**  ISIZE  -  FUNCTION  TO  ESTIMATE  SIZE  OF  STRS  ARRAY 

NO  =  NPFREQ 
IF(NQ.EQ.O)  NQ  =  I 
LADD  =  (MXSTEP-(-NQ-l)/NQ-t-2 
LADDS  =  LSTRS-l-(-4*LADD-f-2*(NPRU  +  NPRS) 

LQUES  =  LADDS  -(-  1SIZE(Q(LE),Q(LMAP),Q(LINDXE).Q(LICTYP). 

Q(LIX),Q(LINDXS)) 

MORMEM  =  LQUES  -  MAXQ 
CALL  DYNAMEM(MORMEM,MAXQ) 

ASSMBLE  LUMPED  MASS  MATRIX  AND  SET  UP  INDEX  TO  STRS  ARRAY  *  * 
LADDS-INPUT  TO  ASS8LE  IS  THE  TOTAL  AMOUNT  OF  STORAGE  IN  Q-ARRAY 
EXCLUSIVE  OF  STRS  ARRAY 

LADDS-OUTPUT  FROM  ASSBLE  IS  THE  SIZE  OF  THE  STRS  ARRAY 

CALL  ASSBLE  (Q(LXC).Q(LYC).Q(LIX),Q(LE).Q{LSTRS),Q(LINDXS), 
Q(LINDXE).Q(LSMASS),LADOS,Q(LMAP),Q(LICTYP)  ) 

WRITE  (6.126) 

•***  LT,LUU,LIXX.LIYY,LTOTAL  ARE  MODIFIED  IN  RESTART  IF  CURRENT 
»♦**  RUN  IS  A  RESTART  RUN 

LPSU  =  LSTRS-l-  LADDS 
LT  .  =  LPSU  +  LADD 
LUU  =  LT  -t-  LADD 
LIXX  =  LUU  -h  NPRU  NPRS 
LIYY  =  LIXX  +  LADD 
LTOTAL  =  LIYY  +  LADD 

►***  PRINT  L-POINTERS 

CALL  PRINTL 

C  ADD  MORE  SPACE  TO  THE  END  OF  BLANK  COMMON.  AT  LEAST  100  WORDS 


C  MORMEM  =  LTOTAL  -  MAXQ 
C  CALL  DYNAMEM(MORMEM,MAXQ) 


onon  nnnnri 
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C 

C****  ZERO  OUT  TIME  DEPENDENT  ARRAYS  XI,  Vl,  Al,  FORCD,  XO 
C 

CALL  ZER0(Q,LX1,LDTN0D-LX1) 

C 

C****  DUMP  Q  AT  THIS  POINT  IF  IQARF  LT.  0  ,FOR  DEBUGGING  RUN. 

C 

IF(IQARF.LT.O)  CALL  DEBUG 
C 

C  IF  THE  USER  HAS  REQUESTED  OUTPUT  TO  THE  COMMON  DATA  BASE  (ITllF 
C  IS  NuN^ERO)  THEN  WRITE  ALL  DISCRIPTIVE  INFORMATION  TO  TAPEll 

IF  (ITllF  .NE.  0)  THEN 

CALL  INIT11(Q(LIX).Q(LXC),Q(LMAP),Q(LIPLT)) 

END  IF 

***•  COMPUTE  TIME  STEP  AND  NUMBER  OF  TIME  STEPS  ***** 

•**•  CALL  FRCIN  W.ZERO  DISPL.  AND  FORCES  ONLY  TO  GET  DELTMX 

CALL  FRCIN  (Q(LXC),Q(LYC),Q(LE),Q{LXl),Q(LVl),Q(LFINT),Q(LSTRS), 
Q(LINDXS),Q(LINDXE),Q(LMAP),Q(LIX).Q(LICTYP).Q(LDTNOD)) 

IF  (DELT  .LT.  0.0)  THEN 
DELT  =  -DELT  •  DELTMX 
ELSEIF  (DELT  .EQ.  0.0)  THEN 
DELT  =  0.5  *  DELTMX 
ENO/F 

IF(TIMEMX.EQ.  0.0)  TIMEMX  =  DELT  *  MXSTEP 
IF(MXSTEP.EQ.  0)  MXSTEP  =  TIMEMX/DELT 
MXC  =  TIMEMX/DELT  -f-  0.001 
MXSTEP  =  MIN0(MXC,  MXSTEP) 

WRITE  (6,211)  DELTMX.DELT.DTMIN.DTMAX, MXSTEP, MXSUB 
IF(DELT  .GT.  OELTMX*0.66667)  WRITE(6,214) 

IF(MITSF.GE.l)  WRITE  (6,212)  (I.Q(LOTGRP-(-l-l),l=l,NUMMAT) 

WRITE  (6,213) 

RETURN 


pormats************* 


1000  FORMAT(  OMITTED  ) 

1002  FORMAT  (//,5X,'CODE:  ',A,5X,’  VRSION;  ’,A,5X,'RUN  ON:  ’,  A, 
5X,’  AT:  ’,  A.  //) 

109  FORMAT  (/////,’  USER  COMMENTS:’,/) 

111  FORMAT  (15X,A) 

110  FORMAT  (A) 

120  FORMAT  (’!’,/,’  TITLE:  ’,A,///) 

130  FORMAT  (7I5.5X,3E10.4) 

5  FORMAT  (’  NO.  OF  D.  0.  F.  (NDGREE)’  ,110/ 

’  NUMBER  OF  NODES  (NUMNP)  ’  ,110/ 

’  NUMBER  OF  ELEMENTS  (NUMEL)  ’  ,110/ 

’  NUMBER  OF  MATERIALS  (NUMMAT)’  ,110/ 
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'  (MO. OF  BOUIMDARY  (MOOES  (NU(MDIS)'  ,110/ 

'  NO.  OF  LOAD  LINES  (NPRES)  ’  .110/ 

’  SLID(NG  INTERFACES  (NSLID)  '  .110/ 

'  MAX.NO.OF  TIME  STEPS  (MXSTEP)'  .110/ 

’  TIME  INCREMENT  (DELT)  ’  .G14.5.1X.A/ 

■  MAXIMUM  TIME  OF  RUN  (TIMEMX)’  .G14.5,1X.A/ 
’  MASS  DAMPING  FACTOR  (ClDAMP)’  .1X.E12.3) 

6  FORMAT  (815) 

126  FORMAT  (///) 

211  FORMAT  (///.5X.  'TIME  STEP  INFO;'./. 


1 

X 

5X,  ' 

MAX.TIME  STEP  COMPUTED  BY  CODE  -  ' 

5X,  ' 

TIME  STEP  USED  =  '.GlO.4,/, 

5X, 

'  DTMIN  =  ',G10.4,/, 

+ 

5X, 

'  OTMAX  =  ',G10.4,/, 

+ 

5X, 

'  TOTAL  NUMBER  OF  TIME  STEPS  =  '  ,110,/ 

+ 

5X, 

'  NUMBER  OF  SUBCYCLES  =  '  .110,/) 

212  FORMAT  (  /.5X.  'MULTIPLE  INTEGRATION  TIMESTEPS  WILL  BE  USED:',/ 

,(7X.I5.G12.4/)) 

213  FORMAT  (//,'  END  OF  INPUT  PROCESSING',/, 

1X.79('*'),/,'1') 

214  FORMAT  (///, 

+  ’  ’./ 

+  '  WARNING,  POSSIBLE  INSTABLILITY.  DELT  SHOULD  BE',/ 

+  '  LESS  THAN  2/3  OF  DELTMX.  OELTMX  IS  BASED  ON  WAVE  SPEED  ',/ 

+  '  CALCULATED  FROM  BULK  MODULUS,  ACTUAL  WAVE  SPEED  IS  PROBABLY  ',/ 
4-  '  HIGHER.',/ 

+  ■ - ’.//) 

END 


n  n 
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SUBROUTINE  READPV(IVOL.THICK.PT,KPRES.LADD,  E.INDXE.lX.BIMP  ) 
C 

c**** 

C****  THIS  ROUTINE  READS  THE  LOAD  LINE  DATA  INTO  ARRAYS  IVOL.PT 
C’***  AND  KPRESS  FOR  USE  IN  ROUTINE  FREEFD 

C***.  IVOL(I)  =  TYPE  OF  LOAD  LINE, 

C****  -2.  AXISYM.PROB.PRESS.LINE 

C****  -1.  PLANE  PROB. PRESS. LINE 

C****  0,  INITIAL  VELOCITY  LINE 

C**’*  1,  FORCE  LINE 

C»***  2.  PRESCRIBED  DISPLACEMENT  LINE 

C****  3.  AXISYM.  NONREFLECTING  BOUNDARY  LINE 

C«**«  4.  PLANE  NONREFLECTING  BOUNDARY  LINE 

C****  |V0L(2)  =  INITIAL  VEL.DIRECTION  CODE, 

C«**«  0,  VEL.IS  NORMAL  AND  TANGENT  TO  LOAD  LINE. 

C****  1,  VEL.IS  IN  GLOBAL  X  Y  DIRECTIONS 

C****  IVOL(3)  =  COMPONENT  DIRECTION  CODE 
C*’**  1,  X-COMPONENT 

C»***  2,  Y-COMPONENT 

C****  IVOL(4)  =  NO. OF  NODES  FOR  CURRENT  LOAD  LINE 

C***.  1V0L(5)  =  NO. OF  TIME,PRESS.PAIRS  {OR  TIME.FORCE  OR  TIME.DISPL.) 

C****  THICK  =  THICKNESS  FOR  PLANAR  PROB.  PRESSURE  LINE 

C****  =  1.0  FOR  AXISYM. PROB.LINE 

C**.*  MPT  =  NUMBER  OF(T|ME.PRESSURE)DATA  PAIRS 

C*»**  PT  =  (T(ME,PRESSURE)PAIR  DATA  **  Tl,Pl,T2,P2,„TNPT,PNPT 

C*.»*  nDNOD  =  NUMBER  OF  NODES  ON  LOAD  LINE 

C****  KPRES  =  LIST  OF  NODES  ON  LOAD  LINE 

C****  ipRES  =  LOAD  LINE  NUMBER 

C*.**  i^aDO  =  STORAGE  USED  IN  STORING  LOAD  LINE  INFORMATION 
«*** 


COMMON/NUM/  NUMNP,NUMEL,NUMMAT,NUMDIS,NDGREE,NPRES,NSLID,MEQ, 
•  DELTMX,LE2,NPLOT 

COMMON/KNTRL/ISYMF.MITSF,IUN5F,ITllF,IMSHF,IHISF, 

+  IQARF.IRSTF 

DIMENSION  KONTRL(IO) 

EQUIVALENCE  (KONTRL(l),ISYMF) 

C 

COMMON/UNITS/LTNM,LLNM,LFNM,LWNM,GMAG,LPNM,LVNM.LDNM 
COMMON/UNITC/  TIMENM,LENGNM,FORCNM,WGTNM,PRESNM  ,VELNM, 

-t-  DENSNM 

CHARACTER*4  TIMENM,LENGNM,FORCNM,WGTNM,PRESNM*12,VELNM*9, 

+  DENSNM*20 

DIMENSION  IVOL(l).  PT(1),  KPRES(l),VOL(l),  E(1),INDXE(1). 

+  IX(10,1),  BIMP(l) 

C 

C****  READ  AND  WRITE  THE  LOAD  LINE  DESCRIPTION  DATA**** 

C 

READ  (5,8)  IPRES. IVOL(4),(IVOL(l),l  =  l,3),INTl,IVOL(5),THICK 
IF(THICK.LE.O.O)  THICK  =  1.0 


n  n  n 
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WRITE(6,9)  IPRES,IVOL(4),IVOL(l),IVOL(2) 

WRITE(6,3)  IVOL(3),INTl,IVOL(5),THICK 
NDIMOD  =  IVOL(4) 

C 

C****  READ  IN  THE  INITIAL  VELOCITY  LINE  DATA  IF  IVOL(l)  =  0**** 
C****  ELSE.  READ  IN  FUNCTION-TIME  PAIRS  **** 

C 

IF  (IVOL(l)  .EQ.  0)  THEN 
READ  (5.10)(PT(J),J  =  1,4) 

WRITE(6,111) 

C 

IF{IVOL(2).EQ.O)  THEN 

WRITE(6.14)  PT(1),VELNM,PT(2),VELNM,(PT(J),J=3.4) 
ELSEIF(IVOL(2).EQ.l)  THEN 
WRITE(6,15)  PT(1),VELNM,PT(2).VELNM,(PT(J),J=3,4) 

ENDIF 
NPT  =  2 
IVOL(5)  =  NPT 
C 

ELSEIF  ((IVOL(l)  .EQ.  J)  .OR.  (IVOL(l)  .EQ.  4))  THEN 
NPT  =  0 
IVOL(5)  =  NPT 

C 

ELSE 

NPT  =  IVOL(5) 

12  =  2*NPT 

READ(5.10)(PT(J).J  =  1,12) 

WRITE(6,100) 

C 

IF  (IVOL(l)  .LT.  0)  THEN 
YVRITE(6,11)TIMENM,PRESNM 
ELSEIF  (IVOL(l)  .EQ.  1)  THEN 
WRITE(6,12)  TIMENM.FORCNM 
ELSEIF  (IVOL(l)  .EQ.  2)  THEN 
WRITE(6,20)  TIMENM.LENGNM 
ENDIF 
C 

DO  1  K  =  l.NPT 
J1  =  2*K-1 
J2  =  Jl-Fl 

WRITE  (6,13)  K,(PT(J),J  =  Jl,J2) 

1  CONTINUE 

*•••  CHECK  TO  SEE  THAT  TIME  POINTS  ARE  IN  RIGHT  ORDER  **•* 

IF  (NPT.GE.2)  THEN 
NN  =  NPT*2-3 
DO  2  J  =  1,NN,2 
IF  (PT(J).GE.PT(J-(-2))  THEN 
WRITE(6,19) 

STOP  'READPVr 
ENDIF 
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2  CONTINUE 
ENDIF 
C 

ENDIF 

COMPUTE  TOTAL  STORAGE  USED  FOR  LOAD  LINE  INFORMATION  **** 

C 

LADD  =  NDNOD  +  2*NPT+6 
C 

C****  READ  IN  NODES  ON  LOAD  LINE 

C****  INDEXING  ON  KPRES  STARTS  AT  II  =  2*NPT  +  1  -  BECAUSE  PT  AND  KPRES 
C****  START  AT  SAME  LOCATION  ***• 

C 

11  =  2*NPT  +  1 

12  =  ll+NDNOD-1 
C 

IF  (INTl.EQ.O)  THEN 

READ  (5,7)  (KPRES(J),J  =  11,12) 

ELSE 

READ  (5,8)  KPRES(ll) 

13  =  12-1 
DO  6  J  =  11,13 

KPRES(J  +  1)  =  KPRES(J)-(-INTl 
6  CONTINUE 
ENDIF 
C 

DO  21  J  =  11,12 

IF(KPRES(J).GT.NUMNP)  THEN 
WRITE(6,23)  KPRES(J) 

STOP 

ENDIF 

21  CONTINUE 
C 

WRITE(6,17)  (KPRES(J),J  =  11,12) 

C 

GET  THE  BOUNDARY  IMPEDENCES  NEAR  THE  NONREFLECTING  BOUNDARY 
C 

IF  (IVOL(l)  .EQ.  3  .OR.  IVOL(l)  .EQ.  4)  THEN 
C 

KK  =  0 

DO  157  K=U,I2-1 
C 

IN1=KPRES(K) 

IN2=KPRES(K-|-1) 

DO  158  l=l.NUMEL 
C 

N=IX(10.I) 

C 

C - 

C  THREE-NODE  TRIANGLE 

C - 

C 


n  n  n  n 


IF  (N  .EQ.  3)  THEN 
C 

DO  159  J  =  l,2 

IF  (IX(J.I)  .EQ.  INI  .AND.  IX(J+1.I)  .EQ.  IN2)  GOTO  700 

159  CONTINUE 

IF  (IX(3,I)  .EQ.  INI  .AND.  IX(l.l)  .EQ.  IN2)  GOTO  700 
GOTO  158 
C 

C - 

C  *••*  FOUR-NODE  QUADRILATERAL  ***• 

C  - 

c 

ELSEIF  (N  .EQ.  4)  THEN 
C 

DO  160  J  =  l,3 

IF  (IX(J,I)  .EQ.  INI  .AND.  IX(J  +  1,I)  .EQ.  IN2)  GOTO  700 

160  CONTINUE 

IF  (IX(4.I)  .EQ.  INI  .AND.  IX(1,I)  .EQ.  IN2)  GOTO  700 
GOTO  158 
C 

C  - 

C  FIVE-NODE  TRIANGLE  •••• 

C  - 

C 

ELSEIF  (N  .EQ.  5)  THEN 
C 

IF  ((IX(1,I)  .EQ.  INI  .AND.  IX(2.I)  .EQ.  IN2)  .OR. 

(IX(2,I)  .EQ.  INI  .AND.  IX(4.I)  .EQ.  IN2)  .OR. 

(IX(4,I)  .EQ.  INI  .AND.  IX(3.I)  .EQ.  IN2)  .OR. 

(IX(3,I)  .EQ.  INI  .AND.  IX(5.i)  .EQ.  IN2)  .OR. 

(IX(5,I)  .EQ.  INI  .AND.  IX(l.l)  .EQ.  iN2))  GOTO  700 
GOTO  158 
C 

c - 

C  SIX-NODE  TRIANGLE  •••* 

C  - 

C 

ELSEIF  (N  .EQ.  6)  THEN 
C 

DO  161  J=1.3 

IF  (IX(J,I)  .EQ.  INI  .AND.  IX(J-|-3.I)  .EQ.  •6.2)  GOTO  700 

161  CONTINUE 
DO  162  3=4,5 

IF  (IX(J,I)  .EQ.  INI  .AND.  IX(J-2,I)  .EQ.  IN2)  GOTO  700 

162  CONTINUE 

IF  (IX(6,I)  .EQ.  INI  .AND.  IX(1,I)  .EQ.  IN2)  GOTO  700 
GOTO  158 


***•  EIGHT-NODE  QUADRILATERAL 


non 
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ELSEIF  (N  .EQ.  8)  THEN 
C 

800  DO  163  J  =  l,4 

IF  (IX(J,I)  .EQ.  INI  .AND.  IX(J+4,I)  .EQ.  IN2)  GOTO  700 

163  CONTINUE 
DO  164  3=5.7 

IF  (IX(J.I)  .EQ.  INI  .AND.  IX(J-3.I)  .EQ.  IN2)  GOTO  700 

164  CONTINUE 

IF  (IX(8,I)  .EQ.  INI  .AND.  ;X(1,I)  .EQ.  IN2)  GOTO  700 
C 

ENDIF 

C 

158  CONTINUE 
C 

700  32  =  INOXE(IX(9.l))  +  5 
DENS  =  E(32) 

YM  =  E(32  +  l) 

PR  =  E(32  +  2) 

BlMP(K  +  NDNOD+KK)  =  SORT(YM*(l.-PR)/(l.  +  PR)/(l.-2.*PR)*DENS) 
BIMP(K  +  NONOO+KK  +  l)  =  SQRT(YM/(2.*(1.+PR))*DENS) 

KK  =  KK  +  1 
C 

157  CONTINUE 
C 

LADD  =  LADO+2*(NONOD  1  < 

ENDIF 

C 

RETURN 


**«*,»**P0RMATS************ 


7  FORMAT  (1615) 

8  FORMAT  (7I5.5X,E10.3) 

9  F0RMAT(//1X,79('*’)/.'  LOAD  LINE  NO.’.IS,/,  '  NO.  OF  NODESM5./, 

’  CONTROL  1  =  M3,’  (TYPE  OF  LOAD  LINE)',/, 

18X,’=  -2,  AXISYMM.PROB.PRESS.LINE’  ,/, 

18X,'=  -1,  PLANE  PROB.PRESS.LINE'  ,/, 

18X,’=  0.  INITIAL  VELOCITY  LINE’  ,/, 

18X,’=  1.  FORCE  LINE’  ,/, 

18X,’=  2.  PRESCRIBED  DISPLACEMENT  LINE’,/, 

1>X,’=  3,  AXISYM.  NONREFL.  BOUNDARY’  ,/, 

18X,’=  4,  PLANE  NONREFL.  BOUNDARY’  ,//, 

.’  CONTROL  2  =  ’,13,’  (NEEDED  ONLY  WHEN  CONTROL  1  =  0  )’,/, 

18X,’=  0,  INITIAL  VEL.IS  NORMAL  AND  TANG.TO  LOAD  LINE’,/, 

18X,’=  1,  INITIAL  VEL.IS  IN  X,Y  DIRECTIONS’,//,) 

3  FORMAT(’  CONTROL  3  =  ’,13,’  (COMPONENT  I.D.  ,  NEEDED  ONLY  WHEN  ’, 
’CONTROL  1  =  1  OR  2  )’,/, 

18X,’=  l(OR  2),  X(OR  Y)  COMPONENT  SPECIFIED’,//, 

NODE  INTERVAL  =  ’,14,’  (FOR  GENERATION  OF  NODES  ON  LOAD  LINE)’,/ 
.’  NO.OF  FN.PTS.  =  ’,14,’  (NO. OF  PIECEWISE  LINEAR  POINTS  USED’, 

’  FOR  PRESS.,’,/, 

23X.  FORCE, ETC.HISTORY  -  NOT  USED  IF  CONTROL  1  =  0)’,/, 
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THICK  =  ’,£10.3.’.  THICKNESS  OF  PLANAR  PROB. PRESS. LINE’,/ 

■) 

C 

111  FORMAT(/  ’  *♦  VELOCITY  SPECIFICATION  **') 

14  FORMAT(/7X,’VEL- TANGENT  ’,1PG20.5,A/ 

7X.’VEL-NORMAL  ',G20.5,A/ 

7X, ’FIRST  NODE  FACTOR’, G20.5/ 

7X,’LAST  NODE  FACTOR  ’.G20.5) 

C 

15  FORMAT(/,6X,'VELOCITY-X  ’,1PG20.5,A/ 

6X,’VEL0CITY-Y  ’,G20.5,A/ 

6X, ’FIRST  NODE  FACTOR’,G20.5/ 

6X,’LAST  NODE  FACTOR  ’,G20.5) 

C 

10  FORMAT  (8E10.0) 

100  FORMAT  (/’  **  FUNCTION  HISTORY  DATA  **’) 

11  FORMAT  (3X,’PT.’,9X,  ’TIME,’  ,A,  2X,’PRESSURE,',  A) 

12  FORMAT  (3X,’PT.’,9X,  ’TIME,’  ,A,  5X, ’FORCE,’, A) 

20  FORMAT  (2X,’PT.’,9X,  ’TIME,’  ,A,  4X,’DISPLACEMENT.’  ,A) 

13  FORMAT  (1X,I4,1PG19.4,G15.4) 

19  FORMAT  (/,5X,’READPV  **  ERROR,  TIME  POINTS  OUT  OF  ORDER’) 

23  FORMAT  (/,5X,’READPV  **  ERROR,  NODE’,15,’  IS  OUT  OF  RANGE’) 

17  FORMAT  (//’  **  NODES  ON  LOAD  LINE  **’  /(1X,20I5)) 

C 

END 
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SUBROUTINE  SOLVE 

.(XC  ,YC  ,XO  ,X1  .VI  ,A1  .SMASS  ,FINT  , FORGO  . 

.  INDXS  ,STRS  .INDXE  ,E  .MAP  .ICTYP  .IX  .OUT  . 

.  NODDIS.ANGLE  .DTNOD  .NPOUT  .PSU  .T  .UU  .IXX  ,IYY  , 

.  ELOUT  ,IPLT  .INDXP  .PR  .INDXL  .SL  ) 

C - 

C  THIS  IS  THE  MAIN  CALCULATION  CONTROLLING  SUBROUTINE.  IT  CONTAINS 
C  THE  TIME  LOOP.  THE  MULTIPLE  TIMESTEP  SUB-LOOP.  AND  THE  CALLS  TO 
C  ALL  THE  MAJOR  CALCULATING  ROUTINES  AS  FOLLOWS; 

C 

C  FREEFD  -  INITIAL  BOUNDARY  CONDITIONS 

C  FREEF2  -  TIME  DEPENDENT  EXTERNAL  LOADS(FORCD) 

C  FRCIN  -  INTERNAL  STRESS  VIA  EOS  ROUTINES 

C  SLIDER  -  SLIDING  INTERFACE  ROUTINE 

C  MOTION  -  EQUATIONS  OF  MOTION  AND  TIME  INTEGRATOR 

C 

C  THE  FOLLOWING  OUTPUT.  POSTPROCESSING  AND  UTILITY  ROUTINES  ARE 
C  ALSO  CALLED: 

C 

C  OUTPUT  -  SELECTED  OUTPUT 

C  ZERO  -  ZEROES  AN  ARRAY 

C  STENGY  -  CALCULATES  ENERGY  TERMS 

C - 

c 

DIMENSION  XC(1)  .YC(1)  .XO(l)  ,X1(1)  .Vl(l)  . 

.  Al(l)  .SMASS(l)  .FINT(l)  .FORCO(l)  ,INDXS(1)  ,STRS(1)  , 

.  INDXE(l)  .E(l)  .MAP(l)  ,ICTYP(1)  .IX(l)  .OUT(l)  . 

.  NODDIS(l),ANGLE(l)  .DTNOO{l)  .NPOUT(l)  ,PSU(1)  , 

.  T(l)  .UU(1)  .IXX(l)  .lYY(l)  .ELOUT(l),IPLT(l)  . 

.  SL(1)  ,PR(1)  ,INDXP(1)  ,INDXL(l) 

INTEGER  PR 
C 

DIMENSION  IQ(1) 

EQUIVALENCE  (Q(l),IQ(l)) 

COMMON  /MAXQ/  MAXQ 
COMMON  /ARRAY/  Q(l) 

C 

COMMON/DYNAM/  DELT.TIME.MXSTEP.NTSTEP.NTSTPl.OMEGA.WMAX.CWAVE. 
-I-  CBULK.CSHEAR.CIOAMP.TIMEMX.DTMIN.DTMAX.MXSUB.POISR 

COMMON/RTIME/  DATE,  RUNTIM,  CREATR,  VRSION 
CHARACTER*8  DATE.  RUNTIM.  CREATR,  VRSION 
C 

C  WARNING - DO  NOT  CHANGE  LOCAT  EXCEPT  AFTER  LTOTAL  OR 

C  SUBROUTINE  DEBUG  WILL  NOT  WORK 

COMMON/LOCAT/  LFINT  ,LXC  ,LYC  .LNODDI.LIX  ,LX1  ,LVl  , 

+  LA1  .LFORCD.LXO  .LDTNOD.LSMASS.LMAP  .LICTYP.LDTGRP.LANGLE. 

+  LINDXE,LE  .LIPLT  .LELOUT.LOUT  .LNPOUT.LINDXP.LPRS  .LINDXL, 

-l-LSLD  .LINDXS.LSTRS  ,LPSU  ,LT  ,LUU  ,LIXX  .LIYY  .LTOTAL 
DIMENSION  LPTR(35) 

EQUIVALENCE  (LPTR, LFINT) 

C 


141 


COMMON/OUTPA/  NPRU,NPRS,NPFREQ,NPIC,NKN,MXN,NPRUS 
COMMON/OUTPC/  TITLE 
CHARACTER  TITLE*80 
C 

COMMON/NUM/  NUMNP.NUMEL.NUMMAT.NUMDIS.NDGREE.NPRES.NSLID.MEQ, 
•  DELTMX,LE2,NPLOT 

COMMON/KNTRL/ISYMF.MITSF,IUN5F,ITllF,lMSHF,IHISF, 

-I-  IQARF.IRSTF 

DIMENSION  KONTRL(IO) 

EQUIVALENCE  (KONTRL(l).ISYMF) 

C 

C 

C****  SET  AND  PRINT  INITIAL  CONDITIONS  WHEN  NTSTPl  =  l*  *  •  *  * 

C 

IF(NTSTP1  .EQ.  1)  THEN 
IF(NPRES.GT.O)  THEN 
DO  244  IP  =  l.NPRES 
LI  =  INDXP(IP) 

NPT  =  PR(Ll+4) 

NNOD  =  PR(Ll  +  3) 

L2  =  LI  +  6 
L3  =  L2  +  2*NPT 

CALL  FREEFD  (XC,YC.X1,V1.F0RCD,STRS.INDXS.IP,NPT, 
NNOD.PR(Ll).PR(L3).PR(L2).PR(Ll+5),  PR(L3)  ) 

244  CONTINUE 
ENDIF 
C 

CALL  OUTPUT(XC,YC.Xl,Vl,Al,IX.E.INDXS,INDXE,OUT,NPOUT.STRS. 
PSO.T,MAP,ICTYP,UU,IXX,IYY,ELOUT,IPLT) 

ENDIF 

C 

CAIL  CPTIME(CPUl) 

C 

C***‘  START  OF  INTEGRATION  LOOP  ♦  *  •  *  • 

C 

DO  260  NTSTEP  =  NTSTPl,  MXSTEP 
C 

C****  MULTIPLE  TIME  STEP  INTEGRATION  LOOP 
C 

DO  1250  NSUB  =  I,  MXSUB 
TIME  =  TIME  +  DTMIN 
C 

C****  COMPUTE  EXTERNAL  LOAD  ♦  •  ♦  •  • 

C 

CALL  ZERO(FORCD,l,MEQ) 

IF(NPRES.GT.O)  THEN 
DO  264  IP  =  l.NPRES 
LI  =  INDXP(IP) 

NPT  =  PR(Ll+4) 

NNOD  =  PR(Ll+3) 

L2  =  LI  +  6 
L3  =  L2  +  2*NPT 
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CALL  FREEF2(XC.YC.X1.V1,F0RCD,STRS,II\JDXS.IP,NPT, 
NNOD,PR(Ll).PR(L3).PR(L2).PR(Ll-r5),  PR(L3)  ) 

264  CONTINUE 

ENDIF 
C 

C****  COMPUTE  INTERNAL  LOADS  VIA  ELEMENT  MODEL  AND  CONSTUTIVE  LAW 
C 

CALL  FRCIN  (XC.YC.E,X1,V1,F!NT,STRS,INDXS,INDXE,MAP, 
lX.ICTYP.DTNOD) 

C 

C’***  SLIDE  LINE  PROCESSING 
C 

IF(NSLID.GT.O)  THEN 
DO  320  JSLID  =  l.NSLID 
LI  =  INOXL(JSLID) 

L2  =  LI  +  11 

CALL  SLIDER(XC,YC,SL(L1),SL(L2).X1,X0.V1,SMASS,A1, 
FINT.FORCD.NODDISJSLID) 

320  CONTINUE 

ENDIF 
C 

C«***  SOLVE  FOR  NEW  DISPLACEMENTS.VELOCITIES  AND  ACCELERATIONS  *  * 

C 

CALL  MOTION(NODDIS,SMASS.XC.YC,XO,Xl,Vl,Al,FORCD.FINT, 

ANGLE, DTNOD) 

C 

12S0  CONTINUE 
C 

C***‘  TIME  IS  RECALCULATED  TO  AVOID  ACCUMMULATION  OF  ROUNDOFF  ERROR 
C****  PARTICULARLY  FOR  PROBLEMS  USING  MITS 
C 

T4M€  =  NTSTEP*DTMAX 

C 

C****  COMPUTE  AND  STORE  ENERGIES  ***** 

C 


CALL  STENGY(INDXE.E,INDXS.STRS,V1,SMASS, ENERGY) 

C 

C****  WRITE  OUT  AND  STORE  RESPONSE  *  DUMP  CORE  IF  REQD.STEP  *  *  * 
C 


C 


CALL  OUTPUT(XC,YC,Xl,Vl,Al,IX,E,lNDXS,INDXE,OUT,NPOUT,STRS, 
PSU,T,MAP,ICTYP,UU,IXX,IYY,ELOUT.IPLT) 


260  CONTINUE 
C 

C****  PRINT  TIMING  SUMMARY 

C 

CALL  CPTIME(CPU2) 

CP  =  CPU2  -  CPUl 

CPI  =  CP/(NUMEL*NTSTEP) 

CP2  =:  CP/(NUMNP*NTSTEP) 

NPl  =  1.0  /  CPI 

NP2  =  1.0  /  CP2 
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WRITE{6,1000)  CPU1,CPU2.NP1.NP2.CP1.CP2 
C 


RETURN 

1000  FORMAT!///, ■  TIMING  SUMMARY:',//. 


-r 

■  SETUP  TIME(S)= 

,1PG12.4./, 

■  END  TIME(S)  = 

G12.4.//, 

4- 

'  ELEMENTS/SEC  = 

16,/. 

-f* 

’  NODES/SEC  =  • 

16.//. 

-f 

’  SEC/ELEMENT  = 

’.  G12.4./, 

-f- 

END 

SEC/  NODE  = 

G12.4) 

Appendix  D 

SAMPLE  INPUT  FILES  FOR  SAMSON2 


l-D  WAVE  PROPAGATION,  DILATATIONAL  WAVE, 


12  6 

1  3 

2  0  200 

0.25 

1  0 

5  0 

0-100 

1 

1  6 

4 

0.00 

1.0 

1.0 

1.0 

0.0  0.00 

1 

0.0 

0.0 

4 

6.0 

0.0 

5 

0.0 

2.0 

8 

6.0 

2.0 

9 

0.0 

4.0 

12 

6.0 

4.0 

1  1 

2  6 

5 

1 

3  3 

4  3 

7 

1 

4  5 

6  10 

9 

1 

6  7 

8  12 

11 

1 

1  21 

-9  21 

4 

1  1 

1  1 

6 

2 

1  3 

2 

1  4  10 

0. 

0. 

2.  8.0E-5 

4.  2.93E-3 

8. 

.01 

10.  1.383E-2 

12.  1.707E-2 

16. 

.02 

35.  .02 

1 

2  3 

4 


SHORT  BAR  W/NRB 


6.  6.17E-3 

14.  1.924E-2 


4 


4 
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1-D  WAVE  PROPAGATION,  SHEAR  WAVE.  SHORT  BAR  W/NRB 


12 

6 

1  12 

2  0 

200 

0.25 

1 

0 

5  0 

0  -1 

0  0 

1 

1 

6 

4 

0.00 

1.0 

1.0 

1.0 

0.0 

0.00 

1 

0.0 

0.0 

4 

6.0 

0.0 

5 

0.0 

2.0 

8 

6.0 

2.0 

9 

0.0 

4.0 

12 

6.0 

4.0 

1 

1 

2  6 

5 

1 

3 

3 

4  8 

7 

1 

4 

5 

6  10 

9 

1 

6 

7 

8  12 

11 

1 

1 

12 

-9 

12 

2 

10 

-4 

10 

6 

10 

-8 

10 

10 

10 

-12 

10 

1111 

6 

2 

1  3  2  2  4  10 

0.  0.  2.  8.0E-5  4.  2.93E-3  6.  6.17E-3 

8.  .01  10.  1.383E-2  12.  1.707E-2  14.  1.924E-2 

16.  .02  35.  .02 

1 

2  3  4  4 

4 
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AXISYMMETRIC  1-D  WAVE  PROPAGATION.  R/L  =  2.0,  SHORT  BAR  W/NRB 


12  6  1  3 

0  0  5  0 

1  6  14 

0.00 

1.0  1.0 

1  58.0 

4  64.0 

5  58.0 

8  64.0 

9  58.0 

12  64.0 

112  6 

3  3  4  8 

4  5  6  10 

6  7  8  12 

1  21 

-9  21 


2  0  200 
0-100 


0.0  0.00 
0.0 
0.0 
2.0 
2.0 
4.0 
4.0 
5 
7 
9 
11 

4 


0.25 


1 

1 

1 

1 


1  1 
6 
2 

1  3 
0. 
8. 
16. 

1 

2  3 
4 


1  1 


0. 

.01 

.02 


3 


1  4  10 

2.  8.0E-5 

10.  1.383E-2 
35.  .02 

4 


4.  2.93E-3  6.  6.17E-3 
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1-D  WAVE  PROPAGATION,  8NQ  ELEMENTS,  SHORT  BAR  W/NRB 
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ANISOTROPIC  1-D  WAVE  PROPAGATION.  DILATATIONAL  WAVE,  SHORT  BAR  W,,  NRB 
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1-D  WAVE  PROPAGATION.  DILATATIONAL  WAVE,  LONG  BAR 
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SMALL  MESH  W/NRB.  SOIL,  SHORTENED  BLAST  LOAD 

49  36  1  7  2  0  4000  O.lOE-5 

00500  -1  000 

1  9  14  YUMA  SOIL  PARAMETERS 


0.00 


0.173E-3 

0.1904ES 

>  0.38 

0.0 

1.0 

0.2644E5 

-0.122 

0.38 

0.4429E5 

-0.156 

0.38 

0.1553E6 

-0.192 

0.38 

0.3442E6 

-0.212 

0.38 

0.2199E7 

-0.999 

0.38 

0.2199E7 

0.0 

0.38 

0.0  0.7249E2 

-0.1E7 

0.65E6 

1 

C.OO 

0.00 

2 

0.00 

6.00 

3 

0.00 

12.00 

4 

0.00 

18.00 
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0.00 

24.00 
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0.00 

30.00 
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0.00 

36.00 
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6.00 

0.00 
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6.00 

6.00 

10 

6.00 

12.00 

11 

6.00 

18.00 

12 

6.00 

24.00 

13 

6.00 

30.00 

14 

6.00 

36.00 
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0.00 
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SMALL  MESH  VV/NRB,  SOIL.  IMPULSE  LOAD 

49  36  1  7  2  0  4000  O.lOE-5 

00500  -1  000 

1  9  14  YUMA  SOIL  PARAMETERS 


0.05 


0.173E-3 

0.1904E5 

0.38 

0.0 

1.0 

0.2644E5 

-0.122 

0.38 

0.4429E5 

-0.156 

0.38 

0.1553E6 

-0.192 

0.38 

0.3442E6 

-0.212 

0.38 
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99999  10  1  1 
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SMALL  MESH  W/NRB,  CONCRETE,  SMALL  BLAST  LOAD 

49  36  1  7  2  0  3000  0.50E-6 

00500  -1  000 

1  9  14  FIBER  REINF.  CONC.  PARAMETERS 

0.05 

0.2516E-3  0.1755E7  0.24  0.01  4.0  1.0  3.0 

0.0  1.0 

0.1125E7  -0.6E-2  0.24 

0.1940E6  -0.32E-1  0.24 
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L^RGE  MESH,  SOU 

.,  SHORTENED  BLAST  LOAD 

256  225 

1  48 

1  0  4000  O.lOE-5 

0  0 

5  0 

0-1  000 

1  9 

14 

YUMA  SOIL  PARAMETERS 

0.00 

3.173E-3 

0.1904E5  0.38  0.01  5.0  l.O  2.0 

0.0 

1.0 

4 

0.2644E5 

-0.122 

0.38 

0.4429E5 

-0.156 

0.38 

0.1553E6 

-0.192 

0.38 

0.3442E6 

-0.212 

0.38 

0.2199E7 

-0.999 

0.38 

0.2199E7 

0.0 

0.38 

0.0  0.7249E2 

-0.1  E7 

0.65E6 
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0.00 

2 

0.00 

6.00 

3 

0.00 

12.00 

4 

0.00 

18.00 

5 

0.00 

24.00 

6 

0.00 

30.00 

7 

0.00 

36.00 

8 

0.00 

42.00 

9 

0.00 

48.00 

10 

0.00 

54.00 

11 

0.00 

60.00 

12 

0.00 

66.00 

13 

0.0C 

72.00 

14 

0.00 

78.00 

15 

0.00 

84.00 

16 

0.00 

90.00 

17 

6.00 

0.00 

18 

6.00 

6.00 

19 

6.00 

12.00 

20 

6.00 

18.00 

21 

6.00 

24.00 

22 

6.00 

30.00 

23 

6.00 

36.00 

24 

6.00 

42.00 

25 

6.00 

48.00 

26 

6.00 

54.00 

27 

6.00 

60.00 

28 

6.00 

66.00 

J. 

29 

6.00 

72.00 

30 

6.00 

78.00 

31 

6.00 

84.00 

32 

6.00 

90.00 

33 

12.00 

0.00 

34 

12.00 

6.00 

35 

12.00 

12.00 

36 

12.00 

18.00 

37 

12.00 

24.00 

38 

12.00 

30.00 

39 

12.00 

36.00 

40 

12.00 

42.00 

41 

12.00 

48.00 

42 

12.00 

54.00 

43 

12.00 

60.00 

44 

12.00 

66.00 

45 

12.00 

72.00 

46 

12.00 

78.00 

47 

12.00 

84.00 
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90.00 

49 
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50 
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53 
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54 
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56 
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42.00 

57 
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48.00 

58 
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59 
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60.00 

60 
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66.00 

61 

18.00 

72.00 

62 

18.00 

78.00 

63 

18.00 

84.00 
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18.00 

90.00 

65 

24.00 

0.00 

66 

24.00 

6.00 

67 

24.00 

12.00 

68 

24.00 

18.00 

69 

24.00 

24.00 

70 

24.00 

30.00 

71 

24.00 

36.00 

72 

24.00 

42.00 

73 

24.00 

48.00 

74 

24.00 

54.00 

75 

24.00 

60.00 

76 

24.00 

66.00 

77 

24.00 

72.00 

78 

24.00 

78.00 

79 

24.00 

84.00 

80 

24.00 

90.00 

81 

30.00 

0.00 

82 

30.00 

6.00 

83 

30.00 

12.00 

84 

30.00 

18.00 

85 

30.00 

24.00 

86 

30.00 

30.00 

87 

30.00 

36.00 
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88 

30.00 

42.00 

89 

30.00 

48.00 

90 

30.00 

54.00 

91 

30.00 

60.00 

92 

30.00 

66.00 

93 

30.00 

72.00 

94 

30.00 

78.00 

95 

30.00 

84.00 

96 

30.00 

90.00 

97 

36.00 

0.00 

98 

36.00 

6.00 

99 

36.00 

12.00 

100 

36.00 

18.00 

101 

36.00 

24.00 

102 

36.00 

30.00 

103 

36.00 

36.00 

104 

36.00 

42.00 

105 

36.00 

48.00 

106 

36.00 

54.00 

107 

36.00 

60.00 

108 

36.00 

66.00 

109 

36.00 

72.00 

no 

36.00 

78.00 

111 

36.00 

84.00 

112 

36.00 

90.00 

113 

42.00 

0.00 

114 

42.00 

6.00 

115 

42.00 

12.00 

116 

42.00 

18.00 

117 

42.00 

24.00 

118 

42.00 

30.00 

119 

42.00 

36.00 

120 

42.00 

42.00 

121 

42.00 

48.00 

122 

42.00 

54.00 

123 

42.00 

60.00 

124 

42.00 

66.00 

125 

42.00 

72.00 

126 

42.00 

78.00 

127 

42.00 

84.00 

128 

42.00 

90.00 

129 

48.00 

0.00 

130 

48.00 

6.00 

131 

48.00 

12.00 

132 

48.00 

18.00 

133 

48.00 

24.00 

134 

48.00 

30.00 

135 

48.00 

36.00 

136 

48.00 

42.00 

137 

48.00 

48.00 

138 

48.00 

54.00 

139 

48.00 

60.00 

140 

48.00 

66.00 

162 

141 

43.00 

72.00 

142 

48.00 

78.00 

143 

48  30 

84.00 

144 

48.00 

90.00 

145 

54.00 

0.00 

146 

54.00 

6.00 

147 
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